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SUMMARY 


An examination of experimental data on the interaction of 
shock waves with a turbulent boundary layer shows that the 
transport of mass and momentum between the free stream and 
the boundary laver are small. Based on this experimental evi 
dence and the characteristics of a turbulent boundary layer, a 
theory is presented to predict the overall boundary-layer changes 
through an interaction. The calculated results compare favor 
ably with the available experimental data 


SYMBOLS 


speed of sound 
parameter defined by Eq. (12a 
average boundary-layer parameter 
momentum flux 
mass flux 
Mach Number 
pressure 
gas constant 
Reynolds Number based on momentum thickness 
temperature 
velocity 
nondimensional velocity 
shock-wave angle 
specific heat ratio 
boundary-layer thickness 
deflection angle of flow 
average boundary-layer parameter 
density 

Subscripts and Superscripts 

= conditions ahead of interaction 

conditions behind interaction 
incompressible 
condition at edge of boundary laver 
condition between two shock waves 
stagnation 
undisturbed 
average condition for boundary layer 
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INTRODUCTION 


i hve BEHAVIOR Of boundary layers in flows with 
pressure gradients has-been the subject of many 
investigations. Some theoretical solutions are avail 
able for both the laminar incompressible and com- 
pressible boundary layers, but there are none for either 
case of the turbulent boundary layer; however, a con 
siderable body of experimental data from measurements 
on the incompressible turbulent boundary layers forms 
the basis for several semiempirical turbulent theories 
Some attempts to extend these theories to compressible 
turbulent layers have been made but the lack of ex 
perimental data hampers these extrapolations. 

An important case in which a supersonic boundary 
layer is subject to an adverse pressure gradient occurs 
when a shock wave approaches a solid boundary. This 
problem has been the subject of considerable experi 
mental investigation. The rapid pressure rise which 
takes place is the distinguishing feature of the shock 
wave-boundary-layer interaction and makes it some 
what different from the usual boundary-layer problem 
The main differences are that the pressure forces acting 
on the boundary layer are much larger than the usual 
values associated with the shear forces at the wall and 
at the outer edge of the boundary layer, and the pres 
sure gradient normal to the wall may be large through 
the interaction region. 

Several analyses have been presented considering 


the problem of separation of a boundary layer on a flat 


surface.'~® Results which correlate reasonably well 
with the interactions in front of forward facing steps 
have been found. However, it has not been possible 
by these techniques to predict the higher pressure ratios 
experienced without separation in the flow in concave 
corners and through incident-reflected shock configura 
tions. 

Because of the large normal pressure gradients, the 


boundary-layer approximations cannot be applied in 
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the interaction region; however, some simplification 
can be made based on the relative magnitude of pres- 
sure forces and shear forces. Some experimental data 
are available which can be used to test the validity of 
these basic assumptions when applied to the various 
interaction configurations. These experimental results 
can be used to determine where the boundary-layer 
assumptions may be applied and the importance of 
shear forces through the interaction region. Based 
on these results, an analytical treatment of these shock- 
wave-boundary-layer configurations can be developed. 


ANALYSIS OF EXPERIMENTS 


The experimental investigations of shock-wave 
boundary-layer interactions have shown several inter- 
esting features. If a boundary layer is caused to sep- 
arate by placing a forward-facing step along a wall, 
there appears to be a characteristic pressure rise asso- 
ciated with this separation. However, if a given pres- 
sure rise is imposed on the boundary layer, by a shock 
wave impinging and reflecting from the wall, pressure 
ratios in excess of the separation pressure ratio found in 
front of a forward facing step can be obtained without 
separation. When separation does occur, it is limited 
to a small region followed by reattachment. Also, for 
the case of an interaction caused by a concave corner, 
pressure ratios without separation in excess of the sepa- 
ration pressure ratio for a forward-facing step seem to 
be possible. Some experimental studies of the flow 
in front of the forward-facing step and through an inci- 
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| 
dent-reflected shock wave are available which 
useful in arriving at an understanding of the mechanjs sas 
involved in these interactions. 

In references 9, 11, 16, and 19, the results of sever 
boundary-layer total head surveys have been present CS 
through some shock-wave - boundary-layer interactioy | 
configurations. These surveys are for forward-facins 
steps at 1/ 2.3, 2.92, and 3.85 and for an incident 
reflected type shock wave interaction at .1/ 2.9) thee 
The first step in interpreting these results is to deter ‘ 
mine the variation of static pressure across the bound std 
ary layer. The static pressure outside of the boundar am 
laver may be inferred from the total head tube r 783 
ings, if the local stagnation pressure is known. Th . 
static pressures at the outside edge of the layer hay |! re 
been calculated under the assumption that the flow y ” 
to the point considered obeved (1) the isentropic | 
and 2 the oblique shock law Che results t suel a 
calculations tor the step and shock-wave case at J/ pbs 
2.9? are shown in Fig. | These results are typical « : 
the other tw \lach Numbers Chis Figure show 
that the static pressure across the boundary laver is : 
approximately constant before and beh nd_ the inter a 
action region but may vary by factors of two or mor sy 
in the interaction region Che static pressure calcu 
lated by considering the stagnation pressure losses 

) bet 
through the shock waves are in much better agreement 
with the wall stauc pressures measured behind the inter oa 


action region, especially in the case of the incident 


reflected shock interaction. 














035 —— 
030 | | | | | | | 
025 sr 
mm | 
a as 020 —_ sf WA 
ae ai i Q Qa 
O15 O~ ; 
hot 
O10 \ 8 
~ Q' 
005 
ot . 
-20 -1.0 6) 
Distance from step face-inches a) 
) 
6 ———_—_— 
030 |- 
>»! 
0285 }-—__—___— i—+ ++ Lp - 
F— , Cie 
ods 020 | / a _| 
oOo — + ++ ar "i 
St OR MOE OS ee” oN 
m, | 
oo} | | | | | = 
005 
Oo Ee 
-10 O LO 
Distance from intersection of incident 
shock and wall- inches b ) 
Fic. 2. Mass flow in boundary layer through shock-wa 
interactions. (a) Forward-facing step, .1/=2.92; (b) incident 


reflected shock wave, / =2.92 


0 deter 
bound 


oundan 


ck-was 


neident 


SHOCK WAVES AND TUBB 
An sis of the mass and momentum in_ the 
boundary layer on each side of the interaction can be 
made based on these boundary-layer surveys. The 
iss in the boundary laver through the interaction 


gion for both the case of the step and shock wave is 


i 
wn in Fig. 2. This Figure shows that the change 
nass in the boundary layer is relatively small. [i 

the momentum in the stream tube comprising the 


iitial boundary laver is calculated, the results shown 


found. The momentum change in the 


Fig. 5 are 











stream tube 1s caused both by shear forces and _ the 
component of the pressure forces acting on the out 
side streamline in the direction along the wall If 
the step type pressure distribution shown in Fig. 5 1s 
issumed, the pressure contribution may be subtracted 
ut These flow models are based on a study of sur 
we pressure distributions and schlieren pictures such 
is those shown in Fig. 4. The corrected curves, in 
Fig. 3 show that there is very litle momentum change 
lue to shear forces in this stream tube. 

[his analysis of the experiments suggests that it 
mav be possible to treat the shock-wave - turbulent 
boundarv-laver interaction in the following way. 
The dissipative flow is considered as a boundary laver 

the regions both before and behind the interaction. 
The equations of conservation of mass and momentum 
between these two regions are satisfied in the absence 

ass or momentum transfer from the free stream 
rthe wall. A flow model of the type shown in Fig. 5 
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is adopted to give the pressure distribution along the 
outer edge of the boundary layer 

This method neglects the skin friction and the mass 
and momentum transfer along the outer edge of the 
boundary layer, but does not restrict the shear forces 
between adjacent parts of the boundary layer. Fur 
ther analysis of the experimental results shows that 
while there appears to be little momentum transfer 1n 
the outer part of the boundary layer, relatively 1m 
portant momentum transfer occurs in the inner part 
of the layer. However, if average parameters are used 
to represent the boundary layer, this shear distribution 


within the layer need not be considered in detail 


BOUNDARY-LAYER REPRESENTATION 


The parameters which will be used to represent the 
boundary layer are those introduced in reference 22 
An average velocity in the boundary layer is defined 
as the ratio of the momentum to the mass # lm 
the 


A quantity «x is now defined as the ratio of ™ to u 


free-stream velocity, « “i u,. A mean density and 
mean temperature may be defined by the continuity 


relation and equation of state 
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the incident shock wave or to the contour of the wal 





The streamline contours are free to adjust themselves 
C) 2) to form a shock pattern or pressure distribution com. } 
patible with the boundary-layer requirements. Foy 
this case, many of the changes do not appear to tak 
place in finite jumps but in a continuous manner 
The interactions spread out over sufficient distang 





7 so that the boundary-layer thickness can change 
For this case, the discontinuity assumption is not ap. } 


(a) 


propriate and no assumption can be made about th 
boundary-layer thickness. 


Many measurements have been made of incompres 
sible turbulent boundary layers under pressure gra- 
dients. These measurements show a tendency for th 


boundary layer to follow a certain class of profile. To | | 
a reasonable approximation, it has been possible t 





describe these profiles by a single parameter. The l 
shear forces within the layer, when given sufficient t 
time, create a boundary-layer profile similar to the 1 
equilibrium profiles. If pressure changes occur very t 


quickly, as in a discontinuous jump, the boundary. | 
layer profile can be radically altered from the equilib- 
rium profile; but when the pressure changes occur 
more gradually, there is a better chance of profiles 











l similar to equilibrium profiles. : 
. ° 1 

t For the case where separation occurs, the boundary 

: TT e 

! Strong layer may have a chance to approach the one-param- 


(c) 


Fic. 5. Shoek-wave boundary-layer models used for cal — 
culating momentum changes in the boundary layer. (a) Forward 
facing step; (b) weak incident-reflected wave; (c) strong incident- 
reflected wave. 
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m = pud = pid/RI iil 

The quantity / is introduced to satisfy the relation Boundary layer 5; | dp 

a ' i. 

T/T, =f — ((y — 1)/2)8@ =f — (v7 — 1)/2]x*w? 

_ ae (a) Normal shock configuration Be 

where w wu a,. The quantities f and « are the two Staats \_* rl 
parameters which will be used to describe the boundary 
layer. 

In order to determine the boundary-layer conditions Boundary layer 
behind the interaction region, an additional relation is = =\— T 
ueeded. The available experimental data give some &, 2, 5, of 
suggestions on what this relation should be. If re 
schlieren pictures of relatively weak interactions are (b) Single oblique shock configuration thé 
examined, Figs. 4(a) and 4(c), 1t can be seen that the 1 bor 
changes appear to take place in finite jumps over most rec 
of the boundary layer, and that the edge of the bound- 
ary layer may be turned, but not displaced, through Be, B, 
these discontinuities. These observations suggest that €2 
these waves may be treated in an analogous way to Streamline 
shock waves in uniform flow if average properties are ® ae 
used to describe the boundary layer. If this type of . 
model is adopted, it is the same as prescribing the 
boundary-layer thickness behind the interaction. Such 5; 
a model might be expected to be correct for the cases Boundary layer} 
of shock-wave interactions of insufficient strength to 1 2 Th 
separate the layer. (c) Incident-reflected 2 shock configuration ! tos 

Once separation occurs, Fig. 4(e), the location of the Fic. 6. Shock-wave boundary-layer models for discontinuity 


shock waves bears no simple relation to the location of analysis 
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SHOCK WAVES AND 


eter family, and this assumption will be made. For 
the weaker shocks, where a discontinuity exists, and 
location of the edge of the boundary layer can be fairly 
well fixed, all empirical assumption of the shape of the 
profile is not necessary. However, for the cases when 
separation occurs, the location of the edge of the bound 


arv layer cannot be specified. 


DISCONTINUITY ANALYSIS 


The discontinuity analysis can be made for both 
normal and oblique shock waves impinging on the 
boundary layer. Equations analogous to the shock- 
wave equations may be developed using the average 
boundary-layer parameters. If the flow 
direction, and velocity at the edge of the boundary 


pressure, 


layer are now assumed equal to the free-stream condi- 
tions, the equations developed can be used to deter- 
mine the boundary-layer parameters f and « behind 


the interaction. 


NORMAL SHOCK 


The general configuration and notation considered 
for the normal shock analysis are shown in Fig. 6(a). 
In terms of the average parameters, the continuity 


equation may be written 


fi - ie = rs2 K1-W, < Ko, 
(1) 
p {fo — | : ae 1) » Ko, 2 KiW, 
a ) { i = * ) ' - 1 ») ~ 
1) 2] | y — 1)/2]e,,2} fx 2/(y 
fir, + 1)/D ley 
] T i — ce 
2y/(y 7 1) 
[) - 
h h Ke) TT 16/17 T 


TURBULENT 
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since 6, 6» for the normal shock configuration 


The momentum equation becomes 


I, + prd ly + pxb» 
and may be reduced to the form 
Pi ae  — ] 2 |K1-% ' 
Xx 

p Uf, — [Cy — 1) 2|woe,,74 

a ¢ 4 a ese ) . 

j Yar U& y l -—j|A2°-U { 2) 

lyKi2w,,.? 4 — [(y — 1)/2]x:2,,29 


If f; and fo, x; and x 1, this relation reduces to the 


classic velocity relation across the normal shock, 


7 + 1)/2]|w,,w, l }) 
By substituting Eq. (4) in Eq. (1) and setting fi, fo, m 
and x: + 
p: > Bes l ? «a! ¥ ] ) 7 + ] 7 
») 
Pi l— i(y — 1)/2]w, 


If the pressure is required to be equal to the free- 
stream value and the boundary-layer properties are 
based on the free-stream velocity, Eqs. (4) and (5) 


may now be used to eliminate w,, and f2 in Eqs. (1) 


and (2), and these two relations may be solved simul- 


taneously for f, and x». When this process is carried 
through, 
‘ 

1)|] (fin ie 

(y—-—1l) (y+ )]i } “73 y= 1)/2Z 
( ) {)) 
f1 — [(y — 1)/2]w,?} 

(l/w ) Ky} Koh )— ] 7) 


hese relations define the boundary-layer parameters behind the normal shock. 


OBLIQUE SHOCK 


The development for the oblique shock is similar to that used in the case of the normal shock. 
of the shock wave within the boundary layer and an average flow deflection angle are defined in the usual way. 


An average angle 


The 


requirement is imposed that the boundary-layer flow deflection angle equals the external flow deflection angle, and 


that the pressure in the boundary layer is equal to the free-stream pressure. 
This general oblique shock configuration is shown in Fig. 6(b). 


boundary-layer shock-wave angle. 


No restriction is placed on the average 
The integration 


required to obtain the average boundary-layer parameters is always performed perpendicular to the streamlines. 


The continuity equation may now be written 
Pp: i. i = 
Pi isc diy o. 29 


and the momentum equation normal to the discontinuity 


Pr V.. BM is Je - hy = 
p» yK1"w,,” sin? Bi + fi — [(y — 


hese equations, with the additional fact that the velocity p: 


to give 


K1}W,, SIN By Kow,, SIN Bo 


[2\(y + 1)]A 


Ki"W ee," Ket’, SIN Pz 

SS) 
Ko-W'..” KyW’,, SIN Oy 
2 |ke"W..” fy — Cy — 1)/2 |i nee 

4) 
2 Kk", fo — I(y — 1)/2 |Ko"e,.° 


irallel to the discontinuity is preserved, may be combined 


= COS-3) 10) 


— [Cy — 1) /(y + 1) |ar2e, 
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Phe unknowns 9), /2, and «2 may be determined by the remaining three equations 
] an ») 
| j 1)/2JWe,- | ; | 
\(y ++ K1"% 2w,.* sin B . 2 
cot 3 ' ; 1) (2 |x,7e,,75 
: | i> 
cot 3 ¥ + le 2% sin- 9 1 1% 2 
; \ kK + 1) /2 |e Ww, sin 
| ) l7 
t yw,, SIN B,w,, Sin w,, sin 
w,, sin 3,w,, sin B - i',.> sin go" 4 | ) cos 
) ) ) Pi 
A 7 4 S11 ) 
v,.7 cos } IS 
kK kK ? j Sil >| 
Eas } | ind 9) provide the relations for 
e external flow and for eiven external shock wave of « agaimst pressure ratio at various Mach Numbers 
Eas () i ind S) mav be used to calculate the his plot shows both the case of single shock and 
ve 1 youndary-laver verage properties. in incident-reflected shock Phe iuitial parameters for 
the reflected shock are taken equal to the parameters 
, behind the incident shocl ind the deflection angles 
CALCULATIONS ‘ 
through the two shocks are taken equal but of opposite 
Using Eqs. (6) and (7) the boundary-laver parameters signs Phis configuration 1s shown in Fig. 6(« The 
behind normal shock wave have been calculated calculation shows that the effect on the boundary laver 
lhe variation of these parameters is shown in Fig. 7 is almost the same for the one or two shock configura 
The calculation has been performed in the region of tions if the same overall pressure ratio is obtained 
lach Number 1.0 to 1.4, the range in which the experi Phis result is similar to the result for the free stream 
ments indicate that this configuration is likely to occur. as might be expected considering the similarity of the 
The initial boundary-layer parameters as a function ot equations, where the difference between one and two 
Mach Number have been chosen equal to the values shock systems is small when the shocks are weak 
lound for flat plate boundary layers shown in Fig. 13. Fig. 11 shows the change in boundary-layer thick 
\ similar calculation has been made for the oblique ness caused by the various shock configurations. A 
shock case, and the results are shown in Figs. S to 11. change in boundary-layer thickness is only possible tor 
n 1] a a e of e poundaryvy-laver parameters is 1@ ODNGCUEe SNOCKS ANd 1S CaAUSCT by l@ Change 1 aire 
he initial value of the 1 1 lay t the oblig hoc] l | by the chat nd 
igain taken from Fig. 13. Fig. S shows the variation tion of the flow Phe boundary laver becomes thinner 
10; and 3, with shock deflection angle @. This Figure in every case in the range considered, except for the 
shows that the boundarv-laver shock-wave angle is normal shock. The calculations have been limited to 
ways somewhat greater than the shock-wave angle this range since the experiments indicate that the 
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TABLE 1 


Properties of the Flat Plate Boundary Layer at Variable 1/7 


M Re, kK Reference 
2.3 17.1 (10%) 0.9276 1.0111 19 
2 92 11.6 (103) 0.9345 1.009 9 
3.85 21.5 (10 0.9393 1.0079 16 
1.97 2.98 (10%) 0.8951 1.018 24, Test 28 
2.6 2.19 (10 0.9187 1.014 24, Test 25 
3.4 4.10 (10 0.9281 1.018 24, Test 20 
1.5 6.59 (10 0.9250 1.011 24, Test 22 


available for com- 
Boundary-layer pro- 


A few experimental results are 
parison with these calculations. 
files are available behind incident-reflected shock wave 
interactions at .\/ 2.92 and total pressure ratios 
of 3.8 2.8. From these experiments, values of 
The experimental 


and 
Ko/K, and 62/6, can be calculated. 
values of these quantities, shown in Figs. 10 and 11, 
compare favorably with the calculated values. 

For the configurations involving large separated 
regions, some means of relating the boundary-layer 
parameters is required. Lacking a sound theoretical 
to fall back on an empirical correlation of these proper- 
ties. A fair body of experimental data exists, especially 
for incompressible flow, that indicates that, to a rea- 
sonable approximation, f and « may be related and the 
boundary layer may be described by a single parameter, 
at least in the regions of mild pressure gradients. Such 
a correlation was first suggested in reference 22, and in 
reference 29 it has been shown that the boundary-layer 
measurements in reference 25 for equilibrium boundary 
layers under three different pressure gradients can be 
described as a one-parameter family with only small 
error. 

This type of assumption has been used in treating the 
problem of shock-induced separation on a flat plate in 
references 1, 5, and 6. References 1 and 6 use differ- 
ent single parameters to describe the separated and un- 
separated layer and reference 5 uses the f-« nomencla- 
ture and a single f-x curve. All of these methods are 
essentially based on the idea that separation can be 
described by a single parameter. Reasonable correla- 
tion was obtained with experiment. 
would lead one to believe that the boundary layer had 
expended itself in overcoming this pressure rise asso- 
ciated with separation, and that it was incapable of 
further pressure rise without influx of momentum. 
This interpretation leaves unexplained the mechanism 


These analyses 


by which the boundary layer is capable of reattaching 
through an additional pressure rise only a short dis- 
tance after separation was occurred, as is observed in 
the incident-reflected type of shock-wave interaction. 
However, an explanation can be obtained of how the 
boundary layer can overcome an additional pressure 
rise without influx of momentum by considering flow 
models appropriate to the configuration being consid- 
ered. In order to demonstrate this point, it will be 
necessary to adopt a one-parameter boundary-layer 
description as was done in previous work. An f-« 
relation, somewhat modified from that used in refer- 
ence 5, which did not fit the boundary-layer profiles 
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measured in shock-wave boundary experiments, will by 
used. 
of Mach Number, it appears to be possible to obtain 
the experimental data 


By assuming the f-« relation to be a function 


a better representation of 
Using these same f-« relations, it will be shown that dif 
ferent flow models can be used to describe both th, 
large separated region on a flat plate and separating 
reattaching flow which takes place in strong incident 
reflected shock-wave interactions. 

The f-« relation which will be adopted for the incom 
pressible flow is that presented in reference 29 based 
23. This curve js 
shown in Fig. 12 labeled JJ = 0. At higher Mach 
Number, there is some data available to show that 
for the zero pressure gradient boundary layer the valu 


on the experiments of reference 


of f decreases and the value of «x increases with Mach 
Number. These data are shown in Fig. 13 along with 
curves for a velocity profile of the form wu, y/6)' 

for values of V = 7 and 9 and listed in Table 1. Not 
enough reliable information is available to draw definite 
conclusions in this matter. For the purposes of the 
present analysis, a set of curves will be constructed 
which are in agreement with the boundary-layer pro 
files measured in the shock-wave boundary-layer ex 
are constructed by trans 


These curves 


forming the incompressible f-« curve by the relations 


periments. 


(1 — x)/(1 — ko) (l — x,)/(1 — k, 


These relations assure that the transformed curves will 
pass through the measured zero pressure gradient point 
1 point, and have a shape similar 


— 


and the« = 1,f = 
to the incompressible curve. 
in Fig. 12 and compared with various boundary-layer 


These curves are shown 


measurements behind shock-wave interactions and in 
the separated region in front of forward-facing steps 
The values used for x and fy are those measured in 
references 9, 16, and 19. These values are used so that 
the results may be compared with the measured re- 
sults of these experiments. These curves are not 
thought to be universal but a reasonable representation 
of the boundary-layer behavior through these shock- 
wave interaction experiments. The comparison with 
the measurements behind the interaction region shows 
that although the argument is not so good as might be 
desired, the present curves represent a reasonable ap- 


proximation of the data. 


SEPARATION ON A STRAIGHT WALL 


To treat the case of separation on a straight wall 
similar to the flow ahead of the forward-facing step, the 
flow model shown in Fig. 5a and found appropriate 
for this case will be used. 

The mass and momentum equations written in the 
direction parallel to the wall in f-« terminology are 
Pokr2W,,62 
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These equations can be combined to eliminate 6, and 6, 


Ko W, 
i — (21) 
kK, W, 


Together with the f-« relation of Fig. 12 and the inviscid 


op) = 
e 


+ ~ 


A 
h Y 


shock relation to prescribe w,, this equation can be 
used to give x2 Kk; as a function of po p;. The maximum 
values of ps p; correspond to the minimum values of 
xk} Therefore, the maximum pressure ratio pos- 


I, + pro; 


yk1-W, - 


6. and 6; may be eliminated to give 


ko[1 — (p p») | 
ky [1 — {p 


‘ 

' 
pi) |) VJ 
For a given incident shock strength, p. and w,, may be 
calculated the shock relation. /p; for 
the configuration shown in Fig. 5(b) should be taken 


from inviscid 
as the value behind the incident shock wave, and for 
the configuration shown in Fig. 5(c) for the value 
associated with the separation pressure rise on a flat 
surface. As stronger shock waves are considered, the 
value p; should be taken equal to the value behind the 
incident shock wave until this value becomes equal to 
the value associated with the separation value. At 
higher incident shock strengths, the pressure p; should 
Eq. (23), 


together with the f-« curve and the inviscid shock rela- 


be kept constant at this separation value. 


tions, is suflicient to determine the conditions behind 
the found, 


62/6, can be found from Eq. (19). 


interaction region. Once f, and x are 


CALCULATED RESULTS 


The separation pressure ratios were calculated using 
Eq. (21) and are shown by the curve in Fig. 14, which 
is compared to the experimental values of references 
20 and The obtained 


through the separation in front of steps has been used, 


28. maximum pressure ratio 
and not the pressure ratio at the point of separation, 
since this analysis should not be applied until after 
conditions have reached their final value and constant 
Static pressure exists across the boundary layer. Since 
it is not known if the initial boundary-layer conditions 
for the experiments of reference 28 are the same as 
those used in the computations, there is some question 
about the comparison with these results. 

Using Eq. (23) and the f-« curve, the results shown in 
been found for the incident-reflected 


Fig. 15 have 
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sible will be the separation pressure ratio. A similar 
analysis was carried out in reference 5; the basic method 
of the present calculation is the same but the numerical 
details are different. 


INCIDENT-REFLECTED SHOCK WAVE WITH SEPARATION 
AND REATTACHMENT 


In the case of the incident-reflected shock wave, a 
similar analysis may be made. Here a flow model of 
the type shown in Figs. 5(b) and 5(c) may be used. 
In both the 
throughout the region between the initial and final 


these models, pressure p; is constant 
pressure rise, so both models are equivalent as far as 
the mass and momentum equations parallel to the wall 
are concerned. 

For this case, the continuity equation is the same as 


Eq. (19) and the momentum equation becomes 


% p YrKoUu 
pb» |} — + 22) 
P = y — 1)/2 |wo" J 
ft 
— 1)/2Ql\e i T YU “ 
99 
yA 
rt e 
~ 1) 2 lw, ie oe “ 


shock-wave case. These results are presented for all 
pressure ratios but should not be correct for pressure 
ratios for which separation and reattachment do not 
occur. Using the f-« curves of Fig. 12, it is found that 
a solution is no longer possible for higher values of inci- 
dent shock strength than those presented. This limit 
occurs at a higher value of x. than the separation value, 
It is 


questionable how much significance should be attached 


which might be expected to be the limiting value. 


to this feature since the /-« curves themselves are ques 
tionable. The important result of this calculation is 
that, by using the same /-« relation, it ts possible to ob- 
tain much higher pressure ratios with the incident- 
reflected shock configuration than the separation pres- 
a flat 
the 


sure ratio associated with the single shock on 


plate. Even though separation occurs within 
interaction region, it is possible for the boundary layer 
to reattach without the addition of momentum from the 
free stream. This result is consistent with the expert- 


mental observations. 


(GENERAL RESULTS AND DISCUSSION 


In the first part of this analysis, the boundary-layer 
properties behind a shock-wave - boundary-layer inter- 
action were calculated by the assumption of a discon- 
tinuity. By this method, both the f and « boundary- 
layer parameters behind the interaction could be pre 
dicted. For interactions involving large separated 
regions, where the boundary layer thickness was given 
time to adjust, an empirical f-« relation based on meas- 


urements of equilibrium boundary layers was used. 
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Both the empirical f-« curves from Fig. 12 and those 
calculated by the discontinuity analysis are shown in 


These curves are of considerably different 


Fig. 16 

shape. The experimental values of f and « behind inci 
dent-reflected shock waves at .1/ 2.92 are shown in 
this Figure. The flow at the lower pressure ratio is 


definitely unseparated and, at the higher pressure 
ratio, separation is just beginning to occur. The point 
for the lower pressure ratio is close to the calculated 
curve while the point for the higher pressure ratio lies 
close to both curves. These two experimental points 
show that the discontinuity relations predict the /-x 
relations up to separation conditions fairly well and 
that the /-« values behind the interactions are consider 
ably different from what might be expected from the 
incompressible equilibrium characteristics. 

Empirical studies of incompressible boundary layers 
suggest that there is a connection between the value of 
x and the separation of the boundary layer. The an- 
alysis presented for the case of separation on a flat 
plate shows that this relation holds for compressible 
boundary lavers if time is given for the boundary layer 
to follow an equilibrium type of /-« relation. It is 
interesting to see whether the same «x values are asso 
ciated with the appearance of separation when the 
boundary-layer parameters follow the calculated /-x 
curves. The values of « calculated by the discontinuity 
analysis for the pressure ratios at which separation is 
first observed experimentally, and the « separation 
values from the empirical curves are shown in Table 2. 
All these values agree surprisingly well, except for the 
last, and indicate that x may be a valid separation cri- 
terion for all turbulent boundary layers. 

In this analysis, two different /-« relations have been 
considered. One relation is based on the equilibrium 
incompressible boundary-layer profiles. This relation 
is for the cases where the changes are gradual and 
the shear forces have sufficient time to reach an equi- 
librium distribution. The other relation describes the 
case behind a discontinuity. This result is obtained 
when the pressure change takes place so rapidly that 
the pressure forces completely dominate the shear forces. 
These two cases represent the extreme cases, and should 
bracket the actual physical changes. 

For the case of a weak shock wave, when the phe- 
nomenon approaches a discontinuity, the discontinuity 
analysis appears to predict the correct /-« curves. For 
the stronger shock waves involving large separated 
regions, the empirical /-« relation is as reasonable an 
assumption as can be made at the present time. The 
important feature of this part of the analysis is that, by 
using the same empirical boundary-layer law, it is 
possible to show that much larger pressure ratios are 
possible with the incident-reflected type of shock-wave 
configuration and an attached boundary laver behind 
the interaction than the pressure ratio associated with 
Separation in front of the forward-facing step. This 
analysis shows the difference in the physical flow model 
which makes the phenomenon possible. 

In this whole analysis, the boundary layer is de 


TABLE 2 
Caleu Calcu 
ed lated 
Experi for for Empirical f 
Configuration mental sep sep Fig. 12 
Normal shock 1.95 0.790 QO S85 0 79 
Incident-reflected 
UV=3 $3 Q S54 S )8 
Incident-reflected 
V/ } 7 ) S55 1.1905 O 8S 
Single ck 
14 1.8 0811 0) 945 0 R80 
Single ck 
VW » § 3.5 0 SIS 1.10] 0.82 
Single shock 
V=33 5.5 75] 1 139 () 89 
* Highest pressure ratio obtainable before separation First 
5, reference last 3, reference 2S 


0 ial 
+ Base on Mach Number behind interaction 


scribed by average parameters. It seems to be satis 
factory to apply the continuity relations to the layer 
as a whole and neglect the shear forces, but this cannot 
be done in detail throughout the boundary layer rhe 
controlling factor in the process is not found in the 
outer boundary layer but in the thin layer near the 
surface. The rate of momentum transport into this 
layer governs the pressure ratio which can be achieved 
before separation and the spreading out of the pressure 
rise in the direction along the wall. As long as only 
average parameters are used, an empirical correlation 
between these parameters and details of the flow near 
the wall will be needed, and the spreading out of the 
phenomenon in the direction along the wall cannot be 
predicted. The next step in a more complete analysis 
would be to include the shear forces in the region near 


the wall. 


CONCLUSIONS 


A theory of shock-wave - boundary-layer interaction 
consistent with available experimental data can be 
formulated on the basis of no external mixing and no 
skin friction. 

For shock waves not strong enough to cause sepa- 
ration, the boundary-layer properties behind the inter 
action can be calculated. 

For strong interactions where large separated regions 
exist, the possibility of the different types of interac 
tions which have been observed can be demonstrated 
by using the appropriate flow models and describing 
the boundary laver as a one-parameter family 

For all configurations, the pressure ratio at which 
separation occurs can be correlated by a boundary-layer 


form factor. 
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SUMMARY tacked in very interesting ways with a direct treatment 
ictior The boundary-laver problem that treats the laminar fre¢ of the partial differential equations by expansion 1n 
rinceton nvective fluid motion produced by a heated vertical circular power series and with the Karman momentum method 
| r which the thert listribution o oO surf: : . k ; 
for which the thermal distri — um the uter surface Consequently, it may be interesting to note that the 
Pres- iries linearly with the distance from the leading edge admits a : i. ; ; 
case where the surface temperature varies linearly with 
2] similarity transformation that reduces the momentum and en ; . ; ; ies 
‘a 1s to two ordinary equations. The exact solutions the height of the cylinder, which is also of practical 
Smooth } tl rv differential system are obtained by numerical importance, admits a similarity transformation 
20-7 methods for parametric values of the Grashof and Prandtl Num 
5 \ comparison is given between the exact solutions and the 
et Pro pproximations that are obtained by the Karman momentum (] MATHEMATICAL FORMULATION OF THE PROBLEM 
using a parabolic assumption for the thermal profile : : ; ” 
- In cylindrical coordinates 7, = the partial differential 
iversity INTRODUCTION equations of laminar boundary-layer theory are 
' 
; [ 
Sharb — : ru } + (714 a 0) | 
hock [ IS WELL KNOWN that the analytical problem asso ) 0 ) O- 1.1 
—— ciated with the laminar free-convective heat trans- n.(Ou. OF) + a.(On. OF 
_ fer from the surfaces of a vertical flat plate can be satis- _ _ 
wart : { ? + (| 7 iM ; + 9f ) 
oni factorily treated by boundary-layer theory. In the v|(o 0 Ou: 0 g3(1 / i 
particular cases where the surfaces are held at constant u.(O7 02) + a.(07 dr) = 
ag temperature, or where certain variations in the surfac a[(027 or?) + (1 7)(OT Op 2 
temperatures are specified, it has been shown!’ that a 
similarity transformation reduces the customary par- where 
ifferential equations of boundary-layer ory toa . . 
tial differenti ul = tions 7 , - _ a layer the t 3 u = vertical component of the local fluid velocity 
ste y re ary ¢ ‘re £ equations. One ol , , : ees 
system of two or lin ry liffe renti il equ tions ¢ fi, = radial component of the local fluid velocity 
Sse dinary equations, whic sse ally describes : ca a . : 
these ordinary equati ns _which essenti lly lescri y = kinematic viscosity of the fluid 
> dynamic: s*havior o : d, is nonlinear; the peer 
the dyn umi¢ il behavi r of the fluid : nonlinez : ¢ = acceleration of gravity 
other, which describes ermal distribution, is + : 
| ther, which describes the Yi : listributi ' 3 = coefficient of cubic thermal expansion 
' linear, although 1 st be realized that the equations z : 
ear, although it an e realized t t | , 1’ = local temperature of the fluid 
are closely coupled. > subsequent s« ons to the z . 
ire clos ly couple 1. rhe subse P on t t T, = ambient temperature 
‘quations are obtained by ‘rical methods since oes 
equatio ire obt uned by numerica et a = thermal diffusivity 
ialytical solutions are unknown. - 
For the closely allied problem of the laminar free- Naturally, the term g3(/° — 79) in the momentum equa 
convective heat transfer from the surfaces of a vertical tion represents the buoyancy effect. An idea of the 


circular cylinder which are held at a constant tempera- physical geometry may be acquired by noting the 
ture, a similarity solution, if extant, has not been ob- schematic diagram in Fig. 1. 
tained; however, the problem has recently been at- Denoting the fixed radius of the cylinder by a, it is 
convenient to introduce the following dimensionless 
” variables and parameters: 
Received June 24, 1957 
i \ preliminary report was presented at the Ninth Interna r=?fa, 2=3/a, Uu, = au,/», 
nels a Applied Mechanics, Brussels, Belgium, Sep- a = abil, 0 = ote TT )/y 
— Scientist, Holloman Air Dev elopment Center and o = v/a, where 9 represents tee Coashot tempera 
eee a a ture and o represents the Prandtl Number. Multi- 
The authors would like to express their gratitude to Frank plying Eqs. (1.1), (1.2), and (1.3) by a», a® »*, and 
Niuman, who programed and obtained the numerical values for gBa* v’*, respectively, in dimensionless terms the funda 
the exact solutions on the Bendix Digital Differential Analyzer mental equations become 
# the Wright Air Development Center, and to Miss Martha 
Tongue, who performed desk calculations and graphical plottings O(ru,) Or] + [O(ru-) Oz 0 1.4 
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Fic. 1 Schematic diagram of the heated vertical cylinder 
showing the coordinate system and notations 





Fic. 3(b Exact solutions for the dvnamic and therm 
files for parametric values of the modified Grashof Number 
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distribution on the 


LAMINAR FREE-CONYN 
u,(Ou, Or) 


O-u,/Or-) + (1/r)(Ou,/Or) +7 (1.9 


Ov Of) + Uu,(OV Or 


l 0 | O-’d Or Z | r)(OvU Or (1.6 


t this point in the analysis it is customary to introduce 


stream function W, such that 


Oy Or ru. and Ow Os — ru (1.4 
one now takes 
y =a2f(r) and 8 = zsh(r)/r' (1.8 


series Of algebraic manipulations shows that the con 
tinuity equation is satisfied identically and that the 
remaining two equations are reduced to the following 


rdinarvy equations 


where primes denote differentiation with respect to / 
9) and (1.10) are ordinary Eulerian differential 


rds 


uations which can be reduced further by the trans 


ration exp ¢ to the following forms 
/ 1) } 
2 F’)F' + 0 l.11 
+h 8)’ + (16 lol F)T=0 (1.12 


where primes now denote difterentiation with respect 
and where // H\) and F(t : The 
transioried to the 


conditions whet 


physical boundar 


. 1 } r 
nensionless variables are 


FiO GO FO 0, J1(0 const BE: 
it the surface of the evlinder and 
F'(o Oo 0 1.14 


in infinite distance fron the evlinder. The bound 


rv condition /7(0 requires that the thermal 


const 


suriace of the cylinder must vary 


linearly with the height of the evlinder 


2) NUMERICAL SOLUTIONS OF THE EQUATIONS 


It does not appear possible to obtain analytical solu 


diiterential 1.11) and (1.12 


ions tor the system 


onsequently, 1t appears desirable to resort to numerical 
obtaining numerical 


procedures The difficulty in 


solutions lies in two facts: the nonlinearity of the first 
liferential equation and the two point boundary cond 

tions hese facts mean that one has to integrate the 
lifferential system as an initial value problem with 


the three known initial conditions and determine the 


two additional initial values F’’(O) and /7’(0), such that 


! 0 and /1(« 0. 
\s a convenience in obtaining approximate values 
F (0) and /7’(0), the Karman momentum method is 


ipplied to Eqs. (1.5) and (1.6 If one denotes r — | 














Exact solutions for the dynamic and thermal ps 
files for parametric values of the modified Grashof 
for Prandtl Number 100.0 


Fic. 3(d 


by y, which represents a dimensionless distance from 
the cylinder wall, and 6 a by 6, which represents the 
dimensionless thickness of the boundary laver, one finds 


by integration that 


bias u--(1 + y)day 
vil + 4 \ Oo! O7/1 ) = 


At a similar point in any analogous treatment it is 
customary to assu:ine dyvnanie and thermal profiles and 
to use Eqs 1) and (2.2) in order to get integral cor 
rections to the constants in the assumed profiles. One 
of the simpler assunptions 1s to approximate the dy 
and to take a pat ib 


namic proile with a particular cubic 


thermal profil 


ness as the 


ola with the same thick 


'.©i, 
ul Ca. V 0 | 0 2.23 
v ee | vy 6 2.4 


These assumptions are very rough approximations 
since they fail to satis'y the differential system when it 
is evaluated at the suriace of the evlinder The as 
suinptions could be improved by taking higher-order 


terms; however, since a first-order result is sufficient 


for later improvement on large scale computers, the 
work that is necessary to obtain higher order approxi 
mations does not appear to be worth while 
Introducing the assumptions (2.5) and (2.4) into 
2.1) and a precisely 


Eqs. and proceeding in 


analogous way to the treatment of the vertical plate, 


one finally sees that 


c-6-(S + 36 j5090)07(4 + 6 1?0¢ 2.0 
( 210 a6 7 T 26 2 6 

1?60 \S > 36 ao AY 6} } + Oo 7 + 6 
I5P00(7 + 26 FT 
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TABLE | 
Comparison of the Exact Solutions With Parabolic Approximations 


Exact Solutions Parabolic Approximation } 
(1) (2) (3 (4 (5 6 
H’'(0 ; max (u; max (u 
o H(O) H'(0) — 4 F’(0 fors = 1 2/6 ( fors =] 
H(O 
0.733 50 107.4 1.852 15.60 2.61 1.58 15.0 » Rg) 
LOO 187.3 ye | 25.99 3.41 1.85 24.9 3.98 
500 511.2 2.978 85.36 8.39 2 68 81.1 Of 
1,000 539.1 3.461 142.7 12.0 3.15 135 ay 
5,000 —4, 802 +. 960 472.2 26.8 $60 444 28.6 A 
10,000 — 18,184 5.818 791.5 38.2 5.43 743 $0.5 
50,000 — 223 ,970 8.479 2630 86.2 8.01 2456 QG).8 
100,000 — 600 ,040 10.00 $415.4 122 9.49 $120 129 ; 
l 100 168.4 2.316 24.70 3.40 2.05 25.2 3.64 | 
1,000 206.2 3.794 135.72 11.0 3.49 137 11.6 ; 
10,000 —24.070 6.407 753.01 34.8 6.04 754 37.0 y 
100,000 — 704, 700 11.05 $202.1 111 10.56 $185 117 . 
10 100 —22.54 4.225 15.74 1.54 3.84 18.5 1.43 : 
1,000 —3,164 7.164 87.06 5.00 6.66 102 $53 
10,000 —83 , 830 12.38 185.0 15.8 11.63 563 14.3 
100 , 000 —1,765,000 21.65 2412.6 50.0 20.5 3151 $5.5 ' 
100 100 346.6 7.466 9.259 0.590 6.74 10.6 0.465 43 
1,000 —8,912 12.91 51.26 1.89 11.8 58.5 1.47 ( 
10,000 — 185,900 22.59 288 .0 6.05 20.8 327 1.67 
100 , 000 —3,556,000 39.56 1600 18.4 36.8 1837 14.8 [ ; 
After 6 is determined for corresponding values of o uct of the coefficient of heat transfer and the radius - 
ut 1 


and %, which is usually done by an inverse process, a of the cylinder divided by the thermal conductivity of pee 


bit of reflection shows that the fluid. For the parabolic assumptions one has 
y o tic 
F"(0) = c (2.8) N, = 2/6 2.11) | 
j 
r > . ‘ ‘ [ > exact s i me has ny 
and, from Eq. (1.8) since 3%) = /7(0), that ind for the exact solution one has Pe 
P H’(0) S 
H'(0O) = (4 — 2/6)H(0) (2.9) N= — 4 2.12 in 
FT(0) 
For sufficiently small values of 6 Eq. (2.7) can be ; . 
: ; Once again, a comparison between the approximate 
solved to give 6 : é 
and exact values of the Nusselt Number occurs in 
5 = (90/7)"/439'/4a-"2(4 + 70)! (2.10) Fig. 2. 
° ° plic: 
gi > re 2c VY ‘ ' ‘ ‘ > 7a , } ‘ 
The values of F"(O) and //’(0) that are given by mins | vn 
Eqs. (2.8) and (2.9) are used as first approximations a : ae 
us. : \ : PI . . ‘Schmidt, E., and Beckman, W., Forsch. Ing-Wes., Vol. | pert 
whose accuracies are improved by a systematic trial- p. 391, 1930. See also Finston, M., Free Convection Past a Ver r 
; va ssi Ma ioe pane ee the 
and-error scheme on a large scale digital computer. tical Plate, J. Appl. Math. Phys. (ZAMP), Vol. 7, p. 527, 1956 
aie - . ‘ J ‘i . ‘i ae : : crea 
lable 1 is a comparison of the approximate and exact * Sparrow, E. M., and Gregg, J. L., The Laminar-Free-Convec- } tl 
‘ . . tio Tr for Fro ) yr Sur ota "ertical Circular ol t 
values; one should compare Column 2 with Column _— Heat 1 ansfer From the Oute Surface of a Vertical Circula 
+3 Cylinder, Trans. Am. Soc. Mech. Engrs., Vol. 78, p. 1823, N cons 
, ‘ P y vember, 1956. See also LeFevre, E. J., and Ede, A. J., Paper Ni \ 
Fig. 2 represents a survey of the initial values as a 1-167, Ninth International Congress of Applied Mechanies _— 
function of H7(0). Figs. 3(a), 3(b), 3(c), and 3(d) are Brussels, Belgium, 1956. 
graphs of the exact solutions for the dynamic and Millsaps, Knox, and Pohlhausen, Karl, The Laminar Free 
thermal profiles for parametric values of the Prandtl Coen Se Themen ree te Senpares fe Fortean Lae F pl 
iG eat Weniiies Cylinder, Journal of the Aeronautical Sciences, Readers’ Forun t 
é ashof } ers. Se ee ra j at 
vane —— : : , . . Vol, 23, No. 4, p. 381, April, 1956 . 
> 2: ‘ ve ~ ¢ ‘ ‘ y ors ° “ ~ ~y I 
In engineering applications an important quantity ‘Goldstein, S. (Ed.), Modern Developments in Fluid Dynan 
is the Nusselt Number .V, which is defined as the prod- ics, Vol. II, p. 641; Clarendon Press, Oxford, England, 1938. 
Re 
i i 
- 7 arrorp TORN Y ~ Y a aa! oan © al Ge TOTTI err ' ‘ 
1958 HEAT TRANSFER AND FLUID MECHANICS INSTITUTE 
atl 
— ; kasi sie ‘ : . ee ; 
Unwersity of California at Berkeley, June 19-21, 1958 f Ac 
Aircr 
The 1958 Institute is devoted to fundamental contributions in thermo and fluid dynamics. Emphasis is given to topics the re 
covering more than one specific field of engineering science; particular attention is devoted to those problems in fluid Caci 
mechanics which require consideration of thermal phenomena (aerothermodynamics), chemical reactions (aerothermo | probl 
chemistry ), and electromagnetic fields (magnetohydrodynamics ). aoe 
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A Pilot Analog for Airplane Pitch Control’ 


N. D. DIAMANTIDES* 


Goodyear Aircraft Corporation 


SUMMARY 


ified pilot analog has been developed that accurately 
e motions of a jet interceptor pilot as he controls a 
irplane. The analog was constructed by use of elec 


mputer components, and its design was clarified by 





stulating mathematical representations of the human dynamics 
ivolv 

The log was tested in a control loop that included simulated 
itch disturbances, the main characteristics of an F-89B jet inter 
eptor, and a movable cockpit mock-up equipped with a simpli 
ied control column and display. The test loop permitted oper 
tion of the analog in parallel with human pilots so the control 
iotions of both, when presented with identical stimuli, could be 
recorded for comparison. Control could be switched, without 
the pilot's knowledge, to analog or pilot. The analog performed 


iccurately that experienced jet pilots often did not immediately 
detect its substitution for themselves in the control loop 

The 
ut it was proved to be representative of an ‘“‘average”’ pilot’s dy 
Moreover, 


iNalog 


log required simple adjustments for each test subject, 


namics within the scope of the simple control task. 
the adjustments permitted measurement of human character- 
isties, thus making the analog a more versatile tool for human dy- 
namics research. Its success encourages the development of more 
mplex pilot simulators that would provide accurate data on 
little-known characteristics of the human as a controller and that 
Iso could be used to predict human performance in specific ma- 


hine systems before costly experimental models are built 


(1) INTRODUCTION 


_ LACK of an accurate method of predicting the 
performance of the human controller often com- 
plicates the development of complex, high-speed mod- 
ern machines. The need for a reliable indicator of the 
performance of the pilot of a high-speed airplane, while 
the machine still is in the early stages of design, in- 
creases proportionately with the increasing difficulty 
of the pilot’s role and with the rising costs involved in 
constructing experimental models. 

An accurate analog of an ‘‘average”’ pilot’s dynamics 
could be put to many uses: 


1) In pilot performance analysis—to test the ap- 
plicability of certain linear and nonlinear analytical 
techniques to the pilot-airplane system. 

to close the loop of the 


») 
{ 


In machine design 


1957. Revised and received December 23, 


Received May 31, 
1957 

1 The airplane mock-up and the original version of the pilot 
were designed by R. J. Mead, the project engineer. The 
project was contractually supported by the U.S. Navy, Bureau 
f Aeronautics, under Contract NOa(s)-54-244-c, with Goodyear 
Aircraft Corporation. The advanced version of the simulator is 
the result of a joint effort of the author and Mr. Mead. Dr. A. J 
Cacioppo was the adviser on the psychological approach to the 
ind E. W. McGraw and M. Mazoh were in charge of the 
Structural design of the mock-up 


ill ilog 


problem, 


“Senior Development Engineer, Avionics and Electronics 


Division 


pilot-airplane system while the airplane still is on 
the drawing board and, thus, permit incorporation 
of changes and improvements prior to construction 


(3) In instrument 
pilot’s use of control information and permit proper 


design—to investigate the 


design of input instruments (gyro-horizon, control 
stick mechanism, radar display) and output devices 
(control-surface linkages, airframe dynamics 

(4) In investigating effects of abnormal aerody 
namic and control-stick-feel parameters on pilot per 
formance—to establish levels of such parameters 
conducive to effective pilot control without over 
taxing pilot physical capabilities. 

(5) In assessing pilot adaptation effects—to de 
termine variations in the day-to-day performance 


of individuals and among pilots. 


A number of pilot simulators have been built and 
tried,' but the novelty of the art and the limited knowl 
edge on which they arbitrarily were based hindered the 
Pilot 
with 


reliable representation of physiological functions. 
simulators, however, when used in combination 
accepted statistical methods, probably constitute the 
only effective method of evaluating the mathematical 
formulas postulated to represent human dynamics. 

In contrast to the electronic analog, the describing 
function and other mathematical approximations of 
pilot physiological characteristics may seem attractive 
as a short method of defining pilot dynamics, but they 
are merely approximations of the analog, which, it 
self, is an approximate device. The tendency to trans 
late human control characteristics into transfer func 
or, where nonlinearities are involved, into de 


is fully justified, however, because 


tions 
scribing functions 
these postulations can be used to formulate a short 
hand overall pilot equation, and the pilot equation 
can be used, though in a limited manner, to represent 
the electronic simulator. 

A pilot equation that can be considered representa 
tive of the simplified electronic pilot analog for air 
plane pitch control, developed by Goodyear Aircraft 
Corporation for the U.S. Navy, is formulated here to 
aid in describing the design of the electronic model. 
This pilot analog was designed for the simple longitud 
nal control of a jet interceptor subjected to forced ran- 
dom disturbances. Thus, a single independent vari 
able was involved, and the task required only con- 
tinuous compensatory control. Considerable refine 
ment of the analog would be required to encompass 
more complex tasks. The analog’s accuracy, how- 
ever, in matching pilot hand motions encourages the 
development of more complex pilot analogs to simulate 
the many variables involved in actual flight. 
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(2) Prtor ANALOG DEVELOPMENT 


Synthetic Approach 

The chief purpose of developing a pilot analog is to 
specify what the pilot is dynamically. By use of the 
available control-engineering data, this problem can be 
approached synthetically or analytically. 

The synthetic approach requires scrutiny of pilot per- 
formance data, as recorded by an oscillograph, on the 
basis of servo experience. It must be remembered that 
this approach requires a good knowledge of pilot con- 
trol behavior and is effective only when the experi- 
menter understands the underlying physiological factors 
as well as the mechanisms used to represent them. 

Taking into account what is known or suspected 
about the physiological make-up of the human con- 
troller, the researcher constructs an analog computer 
circuit or other hypothetical model of the human con- 
trol system. The model of the pilot then is tested in 
the applicable aircraft control situation. The differ- 
ences between human and analog performance outputs 
for identical inputs permit an accurate evaluation of 
the analog and yield data that will prove useful in re- 
fining the analog to the extent that, ultimately, the 
pilot will not detect its substitution for himself in the 
control loop. 

The synthetic approach requires no a priori assump- 
tions regarding the linearity of human dynamics; in- 
stead, it encourages the tendency to compare the hu- 
man control system with familiar discontinuous servo- 
mechanisms. For studies of this type, it is expedient 
to view the pilot as a human servo, or as a nonlinear 
and rather complex control mechanism. Because of 
the considerable evidence that the human operator is 
nonlinear and because the analytic approach primarily 
requires linear approximations, the synthetic approach 
was employed in the development of the pilot analog 
for airplane pitch control. 


Human Control Characteristics 


A disturbance to the airplane pitch attitude was the 
stimulus presented to both the real and simulated 
operators. Perception of this stimulus was the pilot’s 
first task and the first to be considered in designing the 
pilot simulator. 

The human operator normally does not react to 
stimuli unless their magnitude is above a certain thresh- 
old level, set up by the physiological limitations of the 
perceptual organs or by his mental attitude toward 
the situation; the latter is influenced by such factors as 
training and biological make-up. Thus, the operator 
does not sense or does not value small disturbances 
enough to attempt a response. 

Once the pitch deviation exceeds the threshold level, 
the perceived information is transmitted to the brain, 
which decides what kind of corrective maneuver should 
be applied to the controls to offset the disturbance. 
The brain sends the appropriate command to the motor 
elements of the body (muscles and supporting skeletal 
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SCIENCES JUNE, 193958 


parts). The signal flow from sensors through the cep 
tral nervous system to the motor inputs require a meas 
urable time, the average length of which is found to hy 
approximately constant in all normal persons. — This js 
called the reaction time. 

Following perception of the stimulus, the human con 
troller apparently performs two essentially linear oper 
ations: (1) a mental computation of the stimulus, or 
the weighing and summing of position, rate, and accel 
eration to achieve a basis for the decision to act, and 


2) an involuntary placing of definite physical limita 
tions on corrective hand motions, which is caused by 
the neuromuscular, or motor, feedback loop. 

Although seldom, if ever, observed without a reaction 
time delay, mental computation does not require any 
calculation by the pilot. The 
and decision functions are similar to thos 


conscious weighing 
summing, 
used by the rope walker in balancing. 

The type of data used in mental computation can be 
identified with the stimulus receptors involved. The 
movement of the airplane with respect to the horizon 
is transmitted by the eyes as a signal proportional to 
the amplitude of the airplane displacement. Rate-of 
motion information is sensed by the ear canals (viscous 
flow of the endolymph fluid) and the eyes (through 
peripheral vision,” a result of the successive stimula- 
tion of the rods of the macula). Finally, linear and 
angular acceleration stimulates the vestibular organs 
of the inner ear and the proprioceptors of body muscles 
tendons, and joints. It is also conceivable that pres 
sure-nerve nets transmit impact sensations. 

In addition to the components of the original pitch 
error signal and its time derivatives—as modified by 


mental computation, reaction time, and_ threshold 


level—the pilot's estimation of the magnitude of the 
total error signal appears to be an important factor 
in his functioning as a controller. In tests of the pilot 
analog, it was observed that the human pilot initiates 
reflexive control motions concurrently with, or even 
slightly before, zero crossings instead of after one reac 
tion-time interval. The resulting hand motions appear 
on oscillographic recordings as a series of unevenly 
spaced square waves of alternate polarity superimposed 
on the prime control signal and initiated when the pitch 
error signal crosses zero. 

To achieve zero error in the steady state, the pilot 
must center the control stick the instant before the 
airplane attitude reaches the level condition. Without 
prediction, a human controller, functioning similarly 
to a discontinuous servo, could not accomplish such 
centering action in the presence of a finite threshold 
level and reaction time. This discontinuous control- 
signal component and the part of the pilot analog that 


embodies it are called the ‘‘bias.”’ 


Selection of Error Input 


Significant in the design of the pilot analog for ait- 
plane pitch control was the selection of an error input 
that would represent the airplane’s motion parameters 
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PILOT ANALOG FOR 


adequately for pilot perception. To simplify the simu- 
lation of the airplane dynamics, the pitch attitude, 
§,,,, was selected as the variable to supply the pilot 
with cues for control. 

A pilot associates a certain control-stick force per g 
with the operating characteristics of his airplane, which 
include speed (1’), altitude, weight, and center of 
sravity. Therefore, the flight-path attitude, y, can be 
used as the primary input variable because the control- 
stick load-factor increment is | g. 

Control-stick force, as incorporated in the airplane 
mock-up used in these studies, was not affected by the 
flight-path attitude. A mechanism capable of convey- 
ing a true control-stick-feel effect was not employed 
because this would have required simulation of the aero- 
dynamic forces normally acting on control surfaces, and 
it was desired to keep the simulation as simple as 


possible. 


Mathematical Representation of Human Functions 


A block diagram of the pilot analog for airplane pitch 
control is shown in Fig. 1. The order of the various 
operations—mental computation, bias, threshold level, 
is consistent with the order of the 


Although the 


and reaction time 
underlying physiological mechanisms. 
placement of the characteristics may be questioned, 
the circuit based on the design served to trick pilots in 
pilot-pilot analog substitution tests. 

The pitch-error signal, 6,,,, and its two time deriv- 
atives, rate and acceleration, are shown in the mental 
computation section of Fig. 1. The output of the 
computation block is 


nis) = G (S)Ovp(s) = ( e3S- + kos + Ry )Owp(S) (1) 


where s is the Laplace operator representing derivatives 
in the time domain and G, is the describing function for 
mental computation. 

The signal from the mental-computation elements 
results in a threshold output, 12. The describing func- 
tion,’ G,(xv;), relating x. to x), is a function of the ampli- 


tude of This relationship can be stated as 


Xo/X, = G(x) (2) 
rhe describing function for the threshold element (see 
Fig. 2) has no component sensitive to frequency vari- 
ations. Its form for pure sinusoidal inputs, disregard- 
ing all input harmonics that reappear in the output,* 


18 
G(x, — (2/m) sin (a,/x) + 

(2a,/mx;) Vx, —a, (3) 
Where a, is the threshold level and 2/4 rad. per sec. is 


the angular slope of the threshold characteristic curve. 
The will the known 


transfer function 


reaction-time element have 


xa(5)/xels) = CAs) = e€ * (4) 


Where 7 is the time delay. 


Because ‘‘bias’’ is instantaneous (reflexive) and con- 
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veys, primarily, error-direction information, it is ap- 
proximated by a step function of alternate polarity, 
depending on the sign of the error signal. Fig. 3 is a 
diagram for derivation of the describing function for 


bias. The describing function for sinusoidal inputs® is 
X4/X%, = kyG,(x1) = lim [(4k,/2x1) X 
a—>O0 
V1 — (@?/x,7)| = 4R,/axm, (5 


The system kinesthesis, or neuromuscular loop, un 
like the mental computation, threshold, reaction time, 
and bias characteristics, requires the participation of 
the pilot’s biomotor system. The elasticity, damping, 
and inertia common to motor systems with an output 
load (hand and arm) with appropriate active elements 
(muscles) are characteristic of this system. A second- 
order transfer function, therefore, was used to approxi 
mate the forward part of the neuromuscular loop: 


, 


x5//xX3 = k/(ms* + ns + 1) (6) 
where x;/ is the neuromuscular output. 

Eq. (6) does not allow for neural feedback. Nor 
mally, a physiological loop is closed around the motor 
system. Open-loop operation indicates a pathological 
phenomenon called ataxia,* a condition in which the 
effector elements, though normal, are not monitored 
by the central nervous system. Continuous monitor- 
ing of the effectors is centered on the pick-off action of 
proprioceptors, which are rudimentary end-organ re 
ceptors of stimuli caused by the contraction of muscle 
fibers. Proprioceptors are distributed throughout the 
muscles and joints to permit internal feedback of data 


such as control-stick-feel pressure. 


The neuromuscular feedback system involves a 
measurable amount of neural-transmission lag (see 
Fig. 4). The transfer function for kinesthesis, then, is 


V5 


Xs is + 1 
- = G,(s) = k, ) @ 
x Ms? + Ns? + Ls + 1 


The proprioceptors sense differences in pressure or 
force, so the error signal, x;, as well as its associated 
signals, x, and x, (see Fig. 1), should have the dimen 
The element representing 


then should 


sions of force or pressure. 
control-stick-deflection sensitivity, ko, 
provide the error signal with the ratio of control-stick 
displacement to control-stick force so that the pilot's 
motor characteristics can be accurately simulated. 
Frequency components that are not linearly related 
to the stimulus pattern have been observed in the 
human response by Tustin’ and others. Generally, 
such ‘‘noise’’ has been attributed to nonlinear opera 
tions, such as the threshold element. The Navy jet 
pilots who participated in this project, however, ex 
hibited steady hand oscillations, or information feed- 
back,* in varying amounts. The oscillation, super 
imposed on the basic discontinuous servo type signal 
representing pilot hand motion, was nearly sinusoidal 
with a frequency near 1.4 cps. As long as the pitch 
disturbance was anticipated by the pilot, it persisted, 
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even during intervals of no disturbance. The oscillation 
stopped only after the disturbing signal was turned off. 

Because the steady hand oscillations were unmistak- 
ably observed in the performance of jet pilots, the 
methods of piloting high-speed airplanes were examined 
to determine the cause. Piloting requires a knowl- 
edge of the airplane's performance especially at high 
speeds; thus, the pilot must familiarize himself with 
the performance characteristics of his airplane with 
respect to the air mass and other environmental condi- 
tions during high-speed operation. In continuously 
attempting to acquire this knowledge, the pilot imparts 
a succession of short, fast impulses to the control stick. 
Although these oscillations do not deflect the airplane 
from its course, they provide the pilot with a kinesthetic 
sense of the status of his control over the airplane. 

The pilot of a high-speed airplane is particularly 
concerned with changes in the g force, which occur when 
the control stick is deflected or the throttle setting 
changed. Apparently, the oscillatory hand motion 
observed in the performance of the jet-interceptor pilots 
who participated in pilot analog tests was an informative 
feedback process used to inform the pilot of the number 
of gs he could safely withstand at any moment, and, 
therefore, of how much control-stick deflection he could 
afford without danger. The participating pilots, at- 
tributed the oscillatory motion to a desire for continuous 
awareness of their control of the aircraft so they would 
be prepared for the next required corrective maneuver. 
A constant amplitude oscillation of frequency fy = 1.4 
cps was added to the kinesthesis output in the pilot 
analog design to represent informative feedback; this 
oscillation only approximates the phenomenon, but it 
was necessary to duplicate the appearance of the pilot’s 
hand motion even though the effect of informative 
feedback on the pitch error signal is relatively small. 
Too little is known about the oscillatory motion for a 
cross modulation with the pitch-error signal. 

Informative feedback, as it is incorporated in the 
analog, may be represented by the following: 


xe(s) = D[(wo + €)/(s? + wy”) | (S) 
where wy = 2nfy = 8.8 rad. per sec. 
c = function of phase angle 


Based on the equations formulated for the mental 
and physiological functions involved, an equation can 
be derived to represent the overall input-output rela- 
tions of the pilot analog. To the extent that the 
analog’s performance duplicates that of the pilot in a 
given situation (sinusoidal input disturbance), this 
equation can be considered the pilot equation for that 
The equation is derived herein to provide 
Its use for 


situation. 
an analytical description of the analog. 
closed-loop analysis of the pilot-aircraft system is lim- 
ited because of the restrictions on the describing func- 
tion method. 

For the pilot analog shown in Fig. 1, the describing 


function is 
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6. /(S) = Ros,(s) [G,G,(s)G,,(s) + k,G, ]O.cp(s) ak 
D[(wo + 0) /(s? + awe")] 9 


where 6., is the control-stick force and where the ex. 
pressions for the functions G are given by Eqs. (| 
The eight parameters—ko, ki, ke, ks, + 


represent the variable quantities in the 


through (S). 
a,, ky, and D 
pilot dynamic system and have different values for 
different individuals. 


Design of Pilot Analog Circuit 


In contrast to the analytical expression for the pilot 
analog, the simulation circuit is much more general jn 
that it is valid for any input signal, 6,,». 
diagram for the pilot analog is shown in Fig. 5. 

The first pilot operation is mental computation of a 
The threshold characteristic 


The circuit 


linear corrective signal. 
is represented by placing a conventional dead-zone cir 
cuit ahead of an absolute-value circuit, which flips a 
high-gain, prebiased amplifier to a state of polarity 
opposite to its normal state. This amplifier and the 
one immediately following it activate the electronic 
switch, S, when the error signal, 6,,,, exceeds the preset 
(threshold) level. The switch consists of four diodes 
in a series-parallel connection. 
and B are under negative and positive voltages, respec- 
tively, the diodes are cut off ; when the polarities ot these 
points are reversed, the diodes conduct, and the output 


As long as points A 


of the mental computation circuit is fed to the reaction- 
time circuit. 

The reaction-time circuit provides the values of G, as 
stated in Eq. (4). This function® can be approximated 
by 


e “2@TII (1 + 2f,7,s5 + 7,757) + 
n=1 
(1 + 2¢,7,5 + 1,2s7)] = exp | — > uta | 10 


where ¢ is an integer, and its value, when applied to the 
simulation, must be finite and small to avoid com- 
plexity. The approximate simulation of G, is accomp- 
lished by use of three amplifiers as shown in the circuit 
diagram (see Fig. 5). Fig. 6 shows the performance 0! 
the reaction-time circuit for r = 0.20 sec. 

Certain factors limited this simulation of a_ pure 
The most obvious of these is the respons¢ 
The instant a vertical 


reaction time. 
of the circuit to a step input. 
wave is applied to the circuit, a high-frequency transient 
oscillation is initiated; it dies out, however, before the 
time-delay interval has elapsed. The error is of littl 
consequence because of its short duration, its rarity 
(airplane attitude disturbances seldom occur as sudden 
steps), and because the low-pass characteristic of the 
airplane system filters out the greater part of the 
oscillation. 

A second error is significant because, by coincidence 
it is of the same order of magnitude as the neuromuscu- 
lar lag. This error appears as a time lag causing 2 
change in the slope of the wave form that passes through 
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the reaction-time circuit. The feedback impedance of 
the last amplifier in the reaction-time circuit, which 
otherwise would have provided a pure resistance, was 
trimmed to permit the use of this error for the simula- 
tion of neuromuscular lag. 

Both errors are caused by the limited number of 
stages used in the reaction-time simulation. This part 
of the analog circuit is being improved through incor- 
poration of a device for simulating pure time delays 
based on an entirely different concept. 

Bias was represented by a high-gain operational am- 
plifier with an adjustable, limited output. Informative 
feedback was provided by a one-amplifier oscillator 
of the parallel-T, phase-shift type with output ampli- 
tude fixed by a diode-limiter circuit. 


(3) ANALOG CONTROL Loop DEVELOPMENT 


General 


The most illuminating test of the pilot analog circuit 
is its operation in parallel with human pilots in a simu- 
lated flight situation. The system arranged for this 
test included a simulated airplane and a simulated 
environment in addition to the pilot and his analog. 
The environment simulation was restricted to a wind 
disturbance that affected only the longitudinal atti- 
tude of the airplane. The simulated pitch disturbance 
was summed with the output of the airplane simulator 
to provide inputs, 6,,,, to a dynamic cockpit mock-up. 
The pilot in the mock-up received realistic data on the 
disturbance through a display and through angular 
movement of the mock-up seat. He controlled the 
mock-up by means of a control column on his platform. 
The pilot's corrective motions provided inputs to the 
airplane simulator, and thus, completed the _pilot- 
airplane loop. The analog, wired into the computer 
also containing the aerodynamics simulation, was 
activated by the same input signals, @,,,, provided to 
the human pilot. <A single-pole, double-throw switch 
allowed instant switching of mock-up control from the 
pilot to the analog, and vice versa. Fig. 7 shows the 


test control loop. 


Airplane Simulator 

The basic aerodynamic conditions presented were a 
pitch rate proportional to the elevator deflection and 
the damping characteristics of the F-S9B jet inter- 
ceptor. The simulation was simplified by omitting 
the long-period phugoid oscillation of the jet fighter 
because only relatively short-term transients character- 
ized the simple task. It was assumed that the airplane 
maintained a 5,000-ft. altitude and a constant forward 
velocity of 560 ft. per sec. The airplane simulator was 
based on the following transfer function :’ 


6, 1.97(0.73s + 1) 
= (11) 


6, s[(0.24)?s* + 0.245 + 1] 


where 6, is the pitch angle caused by the elevator de- 


flection, 6,. 
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Fig. S shows the circuit used to generate pitch-offse; 
signals continuously on a one-to-one time scale. The cir 
cuit is based on the transfer function 


—(R;/CoRiR3) (ReCos + 1) 
8} RgRsCy Cos? + [Ry + R, + (RR; R;) | Cos 4 lj 


As in the simulation of the pilot, resistor and capacitor 
values were selected to obtain the numerical transl, 
tion for the physical system. 

The cockpit mock-up is a hydraulically driven plat 
form on which are mounted a bucket-type seat, a cor 
trol column, and a cathode-ray tube. The platform 
has one degree of freedom, with a maximum angular | 
deflection of 12°. The bucket seat can be oriented 
face the mock-up axis, for roll motion, or perpendicular 
Angular motions of the 


to the axis, for pitch motion. 
mock-up seat are generated by an electronic analog 

computer as functions of pitch or roll angles. Th 

cathode-ray tube provides a visual display of an arti- | 
ficial horizon. The prime mover is a 1,000-psi, 5-hy 
hydraulic pump (Vickers T3-FAI-262), and the servo 
control valve is a high-quality aviation type (Moog 
V-4). 
rating a 
energizes the control valve. 


A conventional d. c. servoamplifier, incorpo 


100-cps oscillation to avoid valve friction 


The frequency response of the dynamic mock-u 
loop with a 200-Ib. weight in the seat, is essentially 
flat to 1.5 eps and is decreased by 3 db. at 2.5 eps for 
the peak-to-peak input of a 5° pitch angle. The over 
all multiloop system—amplifier, valve actuator, and 
platfiorm—has a transfer function comprising tw 
second-degree and two first-degree polynomials in tht 
Laplacian variable s. These are contributed, respec 
tively, by the platform structure, the hydraulic actu 
ator, the control valve, and the amplifier-input filter 
The amplifier-input filter is by far the most significant 
as far as the response is concerned because of its lov 
cut-off point, which provides a safety measure against | 
the possibility of high mock-up accelerations. 

The effective transfer function for the mock-up 


Diep Op = k/(0.0675s + 1) I; 


where the filter time constant, rz, is 0.0675. Th 
transfer function is incorporated in the simulation (se 
Fig. 5) by means of a condenser in the feedback brane! 
of the first amplifier 

An elaborate instrument simulation, which is chat 
acteristic of most operational flight trainers, was ull 
necessary for this application because the only type ‘ 
performance investigated was the continuous control 
of a single variable. The only visual instrument used 
therefore, was a 5-in. cathode-ray tube masked to thi 
size and shape of a typical gyro-horizon indicator. The 
scope trace, at the horizontal position during level 
flight, was deflected vertically to indicate pitch errors 

The control stick used in the mock-up was a 30-i0 
length of hollow tubing fitted with a jet-intercepto 
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Only one degree of freedom was studied 
the control stick could be 


handgrip 


during a 
placed in either of two positions, corresponding to the 


given run, but 
respective seat orientations. 

Control linkages and constant loads were physically 
simulated in the mock-up by standard airplane cables, 
which were attached to the control stick and loaded by 
idjustable inertial, spring, and viscous drag forces. 
A modified viscous damper of the type used in some 
aircraft was employed. The control-column structure 
did not include a control-stick-feel effect adequate for 
simulation of the aerodynamic forces reflected from the 
control surfaces during maneuvering. Consequently, 
the control column dynamics can be expressed by the 


transfer function 


(as + bs + |) (14) 


G(s) = | 


Experiments with the control column yielded the 
following values for coefficients a and } under a standard 


control column adjustment of O.Sg per inch: 


a = 0.015 sec.* and 6b = 0.44 sec. (15) 


Incorporation of the factors represented by Eqs. (14) 
and (15) complete the simulation. 

The control column dynamics simulation must follow 
the commands of the pilot analog in the same manner 
the actual control stick responds to the commands of 
the human pilot for accurate comparisons of the re- 
sponses of the pilot and his analog. The control- 
stick dynamics are shown in Fig. 5 as the last circuit 
of the simulation. 

The entire simulation apparatus, with the exception 
of the computer and the recorder, was placed in an air- 
conditioned lightproof room for isolation from extra- 
A dim 


illumination of the oscilloscope screen and of the room 


neous disturbances and space reference cues. 
reduced eye fatigue. 


Wind Simulator 


The independent input variable, @,, used to disturb 
the system and produce simulated wind gusts, was ob- 
tained by application of principles similar to those em- 
ployed to simulate noise.» A random wave generator, 
with a sample of radioactive material at its heart, 
emitted particles at a rate defined by a Poisson dis- 
tribution.’ A self-quenching counter tube (Geiger- 
Muller) and associated circuits, energized by the par- 
ticles, produced pulses when the particle-energy ex- 
ceeded a preset level. After amplification, the pulses 
triggered a bistable multivibrator, which generated 
square waves of constant amplitude. The square 
Waves, amplified again in the computer and clipped 
symmetrically to assure a zero mean value, provided 
randomly timed zero crossings. 

Because the zero crossings follow a Poisson distribu- 
tion, the associated power spectrum is easily deter- 
mined and can be altered, by filtering, to make the wind 
disturbance more realistic. Mock-up tests with sub- 
jects experienced in both jet and conventional airplane 
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piloting provided parameter values for the filter, which 


included only one time lag. 


(4) Pitot ANALOG TESTS 


Analog Settings 


The analog parameters seem so numerous that their 
isolation and setting, at first glance, appear to be pro 
hibitive tasks. However, certain simple measurements, 
taken under controlled conditions, provide parameter 
values for at least three analog components: threshold 
level, reaction time, and neuromuscular lag. 

Measurements were made with a test subject in the 
dynamic mock-up. The mock-up was designed so that 
the distance between the subject's eves and the face of 
The 


scope deflection gain was adjusted so that the horizon 


the cathode-ray-tube horizon display was 36 in. 


deflection, with respect to the eye, was one-third that 
of mock-up tilt. Thus, when the mock-up was tilted 
6 deg., the horizon deflection, as shown on the scope, 
was 1.25 in. The display adjustment was established 
on the basis of the subjective judgment of several ex- 
perienced pilots. 
Measurements for the 
mined by subjecting the mock-up to randomly spaced 
The points at which 


threshold level were deter 
step inputs of small amplitude. 
the pilots made corrective responses were established, 
and an average of these points was taken. The aver- 
age was used as a starting point from which the thresh 
old level required for each subject could be set. 

The reaction-time parameters also were investigated 
by observing pilot response to unpredictable step input 
disturbances of amplitudes well above the threshold 
These trials revealed a marked variance in reac 
Reaction-time 


level. 
tion times, even for the same person. 
values were found to follow a skewed distribution," 
with the longer tail of the curve extending toward the 
longer time values. The most probable value emerging 
from such distributions and that generally accepted is 
7 = (0.25sec. However, the task involved in this study 
required a complex of stimuli, presented simultaneously, 
all of which incurred the same response. As a general 
rule,” the more intense the stimulus, the shorter the 
reaction time. Because the stimuli complex presented 
to incur pilot pitch-control motions was equivalent to an 
intense single stimulus, the reaction time of efficient jet 
pilots can be figured as low as 7 = 0.20 sec. 

The experiments in kinesthesis (or neuromuscular 
lag) were not rigorous enough to satisfy the meticulous 
analyst but are adequate for the objectives of the 
A step input was presented to the pilot by the 
The platform remained hori- 
zontal and stationary. A direct feedback loop was 
established between the control column output and 
the display. The step amplitude was held well above 
the threshold level so that the threshold effect could be 
neglected, and it was assumed that the mental-compu- 


study. 
visual display only. 


tation transfer function can be represented by a pure 


gain. To permit this approximation, the rate and 
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acceleration associated with the input were held so 
high, by means of the step signal, that they always were 
above the upper limits of pilot perception. Conse- 
quently, the pilot analog control loop comprised only 
the reaction time, exp (—7s); the control column dy- 
namics, G(s); and the kinesthesis, G,,(s); of which the 
first two components were known. 

After several oscillographic recordings were made, 
al average response curve was established. A supple- 
mentary simulation of the loop (see Fig. 9) was set up, 
and the parameters of the kinesthesis block were deter- 
mined by a trial-and-error procedure in combination 
with curve-fitting techniques. The following values 
of the constants given in Eq. (7) were obtained: / = 
0.05, M = 0.00022, N = 0.00672, and L = 0.0983. 


Analog Performance 

The primary purpose of the experimental and theo- 
retical work that led to the development of the pilot 
analog for airplane pitch control was to evaluate jet 
pilots’ control abilities as affected by the dynamic char- 
acteristics of the airplane and the control column. 
Four jet interceptor pilots recently active in fleet duty 
served as test subjects. The study provided an ex- 
cellent opportunity to attempt the simulation of the 
human link in a closed loop system in which all other 
factors already were represented in a dynamic mock-up 
or electronically. The pursuit of these tasks revealed 
answers to several questions regarding a man-machine 
(pilot-aircraft) system for a continuous rate control 
function. 

The pilots performed routine ten-minute flights. 
Pitch error signals were fed to the pilot analog and to the 
dynamic mock-up servo and display (see Fig. 1). The 
observer could record (1) the pitch error, (2) the pilot’s 
hand motions, and (3) the simultaneously generated 
outputs of the pilot analog. 

Analog parameters (mental computation, bias, and 
informative feedback) could be adjusted until differ- 
ences between the two outputs were negligible without 
interfering with the pilot’s control of the loop. Con- 
trol of the mock-up could be switched to the analog 
without the pilot's knowledge (see Fig. 7). When the 
pilot did not detect the substitution of the analog for 
himself in the control loop, the analog was considered 
representative of his dynamics. 

One might suspect that the human pilot’s output 
wave forms would be too inconsistent to permit precise 
judgment on the nature of the required analog adjust- 
ments. Such inconsistency, however, was seldom ob- 
served in the experienced subjects who participated in 
these tests. Although data on their average accuracy, 
as indicated by the root mean square pitch error, showed 
variance, the modes of operation among pilots were 
Only inex- 


similar and were consistent in individuals. 
perienced subjects exhibited occasional control-stick 
reversals that the analog, having no comparable com- 
ponent for unreliability, could not simulate. 

The five oscillographic recordings, each of which en- 
compasses approximately 40 sec. of ‘“‘flight’’ time, are 
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shown in Figs. 10 through 12 as typical of similar re. 
cordings obtained in more than SO hours of ‘‘flight, 
The recordings illustrate the correspondence of outputs 
of three pilots and their analog circuits for standard 
aerodynamic settings [see Eq. (12) ]. 

The response of the light-airplane pilot (Fig. 10) was 
typical of that recorded for pilots without previous ex. 
perience in jet airplane control. Compared to Jet 
Pilot No. 1 (Fig. 11), he showed little or no informative 
feedback. It is uncertain, however, whether he orig; 
nally exhibited the characteristic and gradually aban 
doned it or whether it was lacking completely, for the 
light airplane pilot participated in mock-up tests befor 
the pilot analog was completed. The informative 
feedback habit is superfluous in mock-up “‘flight’’ and 
eventually is abandoned, even by jet pilots. 

The light-airplane pilot also frequently exceeded the 
mock-up pitch range typical of the jet pilots. Evy; 
dence of his poorer performance was confirmed by his 
efficiency recording (root mean square pitch error 
which showed 25 per cent less smoothing of the random 
disturbance than that achieved by the jet pilots. Pilot 
efficiency was checked by use of the modified analog 
circuit shown in Fig. 13. 

Jet Pilot No. 1 (see Fig. 11) exhibited more informa 
tive feedback than the other subjects. He was mor 
active, also, in other respects, and was the only subject 
who fatigued control-column linkage supports and over 
stressed the mechanism. Despite his high rates of con- 
trol-stick motion, the effects of informative feedback 
are scarcely perceptible in this pilot’s mock-up pitch 
record. Pilot No. 1 showed excellent control of the 
mock-up during the run shown in Fig. 11. 

Jet Pilot No. 2 (see Fig. 12) was the least consistent 
of the four Navy pilots. His informative feedback 
level was less pronounced, and his control efficiency 
was higher than that of the light airplane pilot. The 
recording (Fig. 12) was made with standard aerody 
namic parameters, so estimates of the effects of system 
parameter variations could be based on a comparisot! 
of Figs. 12 and 14. In Fig. 14 aerodynamic damping 
was reduced from the normal value of critical damping 
(0.5) to 0.2. The other parameters remained constant 
The pilot-airplane oscillation, a notably slower oscil 
lation than that of informative feedback, is evident 
The same pilot's hand motion for a flight situation in 
which aircraft sensitivity was reduced 50 per cent 1s 
shown in Fig. 15. The average amplitude of the pilots 
control-stick motion increased to compensate for the 
decreased system gain. Instability was noted but not 
as continuously as in the recording shown in Fig. 14. 

The remarkable similarity of actual and simulated 
pilot hand motions, as recorded by the oscillograph 
proves that the pilot analog system is representative 
of an “‘average’’ pilot's dynamics for the simple control 
task involved. 

The effects of human adaptation, on learning, 
early performance can be observed in tests similar t 
those used to check the pilot analog design. Learning 
effects normally are transient as they occur during that 
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TABLE | 


Typical Pilot Analog Parameter Settings 
Adjustments Required for Test 
Subjects 
(Fraction of Full-Scale Setting 


Light 
Analog Potenti- airplane Jet Jet Jet 

Parameter ometer pilot pilot 2. pilot 8. pilot 4 
Position gain, k, P, 0.530 0.150 0.600 0.600 
Rate gain, k» P 0.300 0.600 0.600 — 0.600 
Acceleration gain 

? P, 0.300 0.400 0.600 0.600 
Threshold level, a; P;, Ps 0.750 0.040 0.075 0.075 
Bias, R; P 0.200 0.300 0.500 0.400 
Informative 

feedback, D P, 0.025 0.100 
Deflection 

sensitivity, & P, 0.100 0.080 0.075 0.080 


initial period of somewhat unpredictable exploration 
when the pilot is familiarizing himself with the machine. 
This type of study was possible because analog param- 
eter settings could be adjusted to suit the individual 
and could be changed to account for day-to-day dif- 
ferences in individuals. (See Figs. 10, 11, 12, 14 and 
bo) 

The naturalness of the control effected by the analog, 
when successfully adjusted, is indicated in Fig. 16. 
After Jet Pilot No. 3 had been “‘flying’’ the mock-up 
for several minutes, control was switched to the analog 
without his knowledge. The pilot, completely un- 
aware of the trick, continued to move the control stick 
in response to the pitch stimuli. 

Table 1 shows parameter settings with the corre- 
sponding analog circuit elements (see Fig. 5). Informa- 
tive feedback data are incomplete because the jet 
pilots, although they carried this habit of their normal 
flight experience to “flight,’’ gradually 
abandoned it as they grew more accustomed to test 
conditions. The gradual reduction of the informative 
feedback oscillation was observed in the responses of 
three jet pilots. The fourth pilot showed no tendency 
toward abandoning it even after eighty ten-minute 
The trait may have disappeared, however, if 


mock-up 


trials. 
the tests had been continued. 

An interesting feature regarding the experimental 
extinction of informative feedback was noted in tests 
of the third jet-pilot test subject. Particular care 
was taken to measure the amount of informative feed- 
back present in his control-stick oscillations when the 
characteristic had been observed consistently in earlier 
trials. The data obtained by measuring adjustments 
in the proper potentiometer of the analog for his day- 
to-day performance were analyzed, and his gradual 
abandonment of the informative feedback habit proved 
to be an exponential function of time.' 


(5) CONCLUSIONS 


Although tests of the pilot analog for airplane pitch 
control were limited by the simplicity of simulated 
conditions and by the relatively small sample of test 
subjects, the analog adequately fulfilled the objectives 
for which it was designed. Its accuracy in matching 
pilot hand motions was proved when active Navy jet 


ABRONAUTICAL SCIENCES JUNE, 1958 


pilots failed to detect immediately its substitution for 
themselves in an applicable control loop. 

The analog was designed to account for pilot adapta. 
tion to new situations as well as for individual differ 
ences. Its usefulness as an aid in human dynannics re 
search, as a result of its versatility in these respects 
was demonstrated, at least for the one task involved 
Evidence of the human pilot’s rapid adaptation char 
acteristics was exhibited when jet-pilot test subjects 
quickly accustomed themselves to the dynamics of the 
system. 

Thus, it can be concluded that the human operator 
while performing a continuous rate-control task (such 
as the control of a pitching airplane), is capable of rapid 
self-adaptation to compensate for a new system dynam. 
ics—at least in instances where no other tasks are in 
volved or where system dynamics do not change 
rapidly. 

By simple adjustment of pilot analog potentiometers 
researchers could determine differences in individuals 
as well as among pilots. Psychologists measured dif- 
ferences in experience and training by measuring the 
adjustments required. It was discovered, for instance 
that a pilot's use of information conveyed by the time 
derivatives of an error increases with experience and 
training.'! 

Analog tests also revealed that the control motions 
of jet pilots for the informative feedback process are 
remarkably consistent in individuals and even from 
pilot to pilot. Thus, the premise that a pilot simulator 
equivalent to any given pilot can be developed, even 
though nonlinearities are involved, is justified. More- 
over, such evidence of pilot consistency indicates that 
an analog computer circuit can be evolved to simulate 
the control motions of an “‘average’’ pilot at a given 
level of training. 

It was found that pilot-induced informative feedback 
is a significant part of the output energy of some jet 
interceptor pilots. Previously attributed to noise 
originating from nonlinear operations, oscillations were 
repeatedly detected in jet-pilot hand motions as a signal 
superimposed on the basic signal of pilot control and, 
after considerable investigation, were proved to be a 
pilot’s means of keeping himself informed of the status 
of his control over the airplane. 

The control efficiency of the human operator also can 
be determined, to some extent, by using a modified 
analog circuit. When the airplane system parameters 
do not vary by more than two from a given set of stand- 
ard parameters, the pilot analog circuit can be em- 
ployed to obtain a first approximation to the root mean 
square control efficiency of the test subject. Higher 
order approximations of pilot control efficiency would be 
of doubtful value unless supported by intensive study 
of the self-adaptation and informative-feedback phe- 
nomena. Thus, the development of airplanes with low- 
pass aerodynamic response involves considerable risk 
unless the source and full effects of these characteristics 
are first determined. 


(Continued on page 394) 
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Theoretical Considerations of Flutter at High 
Mach Numbers 


HOMER G. MORGAN,* HARRY L. RUNYAN,** ann VERA HUCKEL*** 
Langley Aeronautical Laboratory, NACA 


SUMMARY 


ne of the theories for two-dimensional oscillatory air forces 


Som 
which may be applied in flutter calculations at high Mach Num 
bers are discussed. These include linear theory, Van Dyke's 
second-order theory, piston theory, Landahl’s method, tangent 
wedge and tangent-cone approximations, Newtonian theory, and 
a new nonlinear-pressure method. A comparison of the theories 


is made by showing the results of flutter calculations for Mach 
Numbers up to 10, and the possibility of flutter at these higher 


Mach Numbers is pointed out 
Results of flutter calculations are shown to illustrate the vari 
ous effects arising from a nonlinear thickness theory The possi 


bility of large flutter speed thickness effects which depend on fre 


quency ratio is shown 
and flutter speed trends with center of gravity and elastic axis 


The influence of airfoil shape is discussed 


locations are presented. Some possible refinements of piston 


discussed for use at very high Mach Numbers. These« 


theory ar« 
include the use of local flow conditions and the use of Newtonian 
theory over the leading edge of a blunt-nosed airfoil 
SYMBOLS 
speed of sound, ft./sec 
semichord, ft 
chord, ft 
C. = lift coefficient 
( = pressure coefficients, (p p pV2/2 
instantaneous airfoil-surface location, ft 
h wing vertical displacement at elastic axis, ft 
ie = mass moment of inertia of wing section per unit 
width, lb.-sec.? 
= reduced frequency, bw/ V 
L:, M = 1, 2, 3, 4) oscillatory total lift and moment 
derivatives defined as in reference 2 
n = mass of wing section per unit width, Ib.-sec ?/ft.* 
V = Mach Number 
= pressure, lb./sq. ft 
t = free-stream pressure, lb./sq. ft 
a = I.,/mb? 
t = time, sec 
V = velocity, ft./sec 
downwash velocity, ft./sec 
Cartesian coordinates, ft 
x = distance from elastic axis to center of gravity, 
positive aft, fraction chord 
\ = elastic axis location, fraction chord 
} = airfoil surface, measured from airfoil mean line, ft 
a = angle of attack, rad 
p = V lf? l 
= ratio of specific heats 
u = mass ratio for a wing section, 1 /4b?p 
im = mass ratio for a cone, (total mass/8)*p) 
p = air density, slugs/ft 


Presented at the Aeroelasticity—II Session, Twenty-Sixth An- 
nual Meeting, IAS, New York, January 27-30, 1958 

* Aeronautical Research Engineer. 

** Head, Dynamic Loads Analysis Section 
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airfoil thickness, fraction chord 


r = 
o = velocity potential 

6 = cone semiapex angle, rad 

w) = circular bending frequency, rad./sec 
War = circular torsion frequency, rad./sec 
@ = 2kM?/p? 


INTRODUCTION 


hn STRUCTURES of high-speed aircraft and missiles 
are, of necessity, designed for maximum load 
carrying efficiency. This implies minimum weight to 
sustain a given load and often results in a structure 
which is relatively flexible, thus inviting such aero 
elastic problems as flutter and divergence. Also, cer 
tain nonlinear effects in the aerodynamics, of com 
paratively small importance at subsonic speeds, can 
have destabilizing influences on these aeroelastic 
problems at high speeds. 

Unsteady aerodynamic theory and the oscillatory aero 
dynamic forces needed for flutter calculations have re 
ceived much attention in the supersonic speed regime 
Correspondingly, little effort has gone toward extending 
this work into the high supersonic and hypersonic speed 
regions, and few applications to flutter have been made. 
However, there are several methods available or adapt 
able for calculating unsteady aerodynamic forces in 
this speed range. 

It is the purpose of this paper to discuss and to at 
tempt to bring into perspective some of the aerodynamic 
theories that are available at high supersonic speeds 
and to show their application to flutter calculations 
The theories to be discussed are the Van Dyke second 
order theory, piston theory, linear theory, the Landahl 
extension of piston theory, Newtonian theory, tangent 
wedge approximation for hypersonic flow, and a new 
theory based on the use of the nonlinear-pressure ex 
pression. 

These various theories will be compared and dis 
cussed with reference to their adequacy for prediction 
of aeroelastic effects. In addition, some flutter trends 
at high Mach Number will be discussed. 


HiGH MacH NUMBER AERODYNAMIC THEORIES 


Before discussing some of the high Mach Number 
aerodynamic theories which are available for aeroelastic 
calculations, some general comments will be made. 

The basis for many supersonic flutter calculations 
has been the work of Garrick and Rubinow.! This 
method has proved to be very useful in the low super- 
sonic speed region and is based on linearized two- 
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dimensional aerodynamic theory. Linearized theory 
has also been used to study three-dimensional aspect 
ratio and plan-form effects in the supersonic region. 
For example, references 2 and 3 treat special plan 
forms while reference 4+ presents a method for handling 
an arbitrary plan form. 

Several attempts have been made for an improved 
supersonic aerodynamic theory which will account for 
all or part of the nonlinearities not treated by linear 
theory. These nonlinearities appear in the aerody- 
namic forces as the result of airfoil shape, thickness, 
boundary layer, Reynolds Number, finite shock waves, 
ete. Although the general problem is largely untouched, 
some of these effects have been studied. For example, 
Carrier's work in reference 5 examines the effects of 
the bow shock wave on an oscillating wedge. Several 
methods are available which include an effect of airfoil 
shape and thickness. For example, Jones and Skan® 
used a numerical integration technique to predict forces 
and moments on a biconvex airfoil. Van Dyke’ found 
a second-order solution for unsteady flow by a method 
of successive approximations. His solution is related 
to linear theory for the oscillating case as Busemann’s 
solution is related to that of Ackeret for the steady 
case. 

To reconcile certain discrepancies which existed be- 
tween the results of Van Dyke and some previous 
writers, Lighthill,’ reasoning from a suggestion by 
Hayes,* presented a check at high Mach Number 
which verified the solution of Van Dyke. This pro- 
posal by Lighthill, since called piston theory, has 
proved to be a very valuable flutter tool in itself and 
has been discussed in this light by Ashley and Zar- 
tarian in reference 10. Piston theory is especially at- 
tractive for flutter calculations because of its simplicity 
in comparison with other supersonic theories. _ Piston 
theory not only agrees with linear theory at high Mach 
Numbers but it is a bona fide high supersonic theory 
which agrees with the more exact second-order solu- 
tion of Van Dyke when thickness terms are included. 

By assuming large Mach Numbers, Landahl'' has 
developed another nonlinear method by use of an iter- 
ation technique similar to that of Van Dyke. Second- 
order piston theory is a special case of this work. 

Still another technique available for accounting for 
a portion of the nonlinearities of the supersonic flow 
problem involves using the linear velocity potential 
but retaining the pressure coefficient in its general 
second-order form. This technique will be referred 
to as the nonlinear-pressure method and will be dis- 
cussed later. 

The methods just discussed are to be applied in the 
supersonic and high supersonic speed regions. For 
purposes of this discussion, 1J ~ 4 can be defined as 
being hypersonic. Piston theory, as originally sug- 
gested by Hayes in reference 8, is a hypersonic theory 
for slender profiles. The tangent-wedge approxima- 
tion, discussed in reference 12, is also a hypersonic 
theory, similar to piston theory, and can be adapted 
for unsteady flow. Each element on the surface is 
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assumed to be subjected to the pressure which would 
act on a semi-infinite wedge whose half-angle is the 


local slope. The pressure relation is derived by con. 


sidering the oblique shock relations for large Mach | 


Numbers and small angles. 

Another procedure which involves the Newtonian 
approximation may also be used at hypersonic speeds 
The flow region is assumed to be an infinitely thin 
layer lying between the shock wave and the airfoil 
Local pressures are determined by the change in fluid 
momentum as it follows the airfoil contour. 

The previously mentioned methods are for two. 
dimensional flow. A second-order three-dimensional 
theory is given by Leslie in reference 13, based on an 
iteration method similar to Van Dyke's. It is applied 
to a rectangular wing with a specialized thickness dis. 
tribution, and static lifts and moments are calculated 
The possibility is presented that this technique could 
be extended to the quasi-steady case for oscillatory 
derivatives. 

In order to understand the relation between the 
various aerodynamic theories, some of those found 
useful for flutter calculations will now be discussed 
in more detail. Linear theory will not be reviewed 
since a large amount of literature is available. 


Van Dyke’s Method 

Van Dyke in reference 7 has presented a solution 
for the oscillating two-dimensional airfoil which is 
correct to second order. This solution contains the 
effects of airfoil shape and thickness on the forces and 
moments acting on the airfoil. By obtaining a solu- 
tion correct only to second order, Van Dyke is able to 
work with the velocity potential equation. (Shock- 
wave effects enter as third-order terms so that to third 
order and higher a potential will not exist.) The 
second-order nonlinear perturbation equations which 
result are, for the velocity potential, 


Gy — B°brr — 2(M?/V) br — (M?/V*)bu = 
M*}(y — 1) [be + (1/V) oe] (bez + by) + 


' 


26-622 + 2b,br, + (2/V) (b:621 + 6,6,0! (1) 
and for the pressure coefficient, 


C, = —26, — (2/V)6, + 6%." — ¢,* + 
2(M?/V)drdt + (M?/V2)be? = (2) 


The boundary conditions require that the flow be 
tangential to the airfoil surface. 

The more familiar linear theory presented by Gar- 
rick and Rubinow, reference 1, results when the right 
side of Eq. (1) is set equal to zero and the product terms 
in the pressure coefficient, Eq. (2), are dropped. In the 
linearized theory, the tangency condition is satisfied 
at the mean line of the airfoil. 

Van Dyke solves Eq. (1) by an iteration technique 
using the solution to the linearized equation as the first 
step. The linear solution is substituted in the right 
side of Eq. (1), giving a linear, nonhomogeneous differ- 
ential equation which can be solved with the correct 


boundary conditions. Pressures are then found at the 
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FLUTTER AT HIGH 


irfoil surface from Eq. (2). Lift and moment coeffi- 


cients are presented in reference 7 expanded in powers 
of a frequency parameter, up to and including the cube 
of frequency. In principle, the solution can be found 
for arbitrary frequencies, and reference 7 gives the re- 
sult for a wedge oscillating about its vertex. How- 
ever, the labor involved in finding the lift and moment 
on a general airfoil at arbitrary frequencies makes 
such a procedure extremely difficult. Nevertheless, 
the coeffcients as presented to the cube of frequency 
should cover most flutter cases of practical interest. 
Inasmuch as this method is based on the general 
potential flow equations, it is expedient to use this 
second-order solution as the standard of comparison 


for other theories at supersonic speeds. 


Nonlinear-Pressure Method 


In Van Dyke's second-order solution, nonlinearities 
in the aerodynamic forces arise from two sources—the 
solution of the nonlinear velocity potential equation 
and from the nonlinear pressure coefficient. If it is 
assumed that contributions to the thickness terms by 
the nonlinear part of the velocity potential are small, 
a thickness effect can be found by using the linear 
velocity potential with the nonlinear pressure coeffi- 
cient. 

The linear velocity potential—1.e., 
Eq. (1) with the right side of the equation set equal to 


the solution to 


zero iS 


vealed [((@/cM) (w — &)] dé (3) 


where w(é, ¢) is the downwash at the mean line of the 
airfoil and the boundary conditions are satisfied as y 
approaches zero. However, downwash w(é, ¢) can be 
interpreted as the downwash at the surface of the 
airfoil without too much error, provided the thickness 
remains small, and will be, for a pitching and translating 


airfoil, 


+ lh + Va+ (i - 2bxo) a] + 


jupper surface 


nan 


VY’(é) s 
(lower surface 
where }(£) is the equation which describes the airfoil 
contour as illustrated in Fig. 1. Using Eq. (4) for the 
downwash, the velocity potential, Eq. (3), can be ex- 
panded in powers of a frequency parameter. This ex- 
pression is used to get the pressure coefficients, using 
nonlinear Eq. (2) and retaining only terms which are 
linear in displacement. The lower 
surface can then be subtracted from the pressure on the 
upper surface and thickness terms retained. 

Such a procedure obviously will not yield the com- 
However, as will be shown 


pressure on the 


plete nonlinear solution. 
subsequently, this method does give the major portion 
of the thickness effects. The primary advantage of 
this technique is that frequency effects are readily ob- 


tained. While Van Dyke’s method is currently avail- 
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able only to cube of frequency, the nonlinear pressure 
approach can be used to obtain an arbitrary number of 
frequency terms without prohibitive labor; while cubic 
frequency terms will prove adequate for a large number 
of flutter problems, other cases, especially in the lower 
supersonic region, will occur where more frequency 


terms are needed. 


Piston Theory 


Piston theory as an aerodynamic tool has been dis 
The 


basic assumption of piston theory is that the local pres 


cussed in reference 10 by Ashley and Zartarian. 


sure on an airfoil is related to the local fluid velocity 
normal to the free stream in the same manner that 
the pressure on a piston in a one-dimensional channel 
is related to the velocity of the piston. For this assump 
tion to be valid, Mach Number must be much greater 
than 1.0. Pressures are time dependent only through 
the boundary conditions—that is, through the normal 
fluid velocities at the airfoil surface. 
point are assumed independent of the conditions at 
In the actual case when J/* > 1.0, 


Pressures at a 


any other point. 
since the component of fluid velocity along the x axis 
changes very little and is much greater than the speed 
of sound, disturbances cause practically no change in 
the x component of velocity from one point to another. 

The piston theory third-order expansion for the pres 
sure at a point, as given in references 9 and 1(, is 


Pp — po = p21 (w/a) + [(y + 1)/4] (w/a)? + 
((y + 1)/12] (w/a)*} (5) 
where w is the local downwash, given by Eq. (4) for a 
pitching and translating airfoil. Lighthill, in refer- 
ence 9, points out that Eq. (5) can be used with good 
accuracy provided \/? > 1 and wa < 1. Since the 
downwash w must be less than the sound speed, a, a 
limitation is placed on amplitude of oscillation and air- 
foil thickness. Since in flutter small perturbations are 
considered, this is primarily a limit on the airfoil thick- 
ness. This also implies that piston theory will not be 
applicable near the leading edge of blunt-nosed airfoils 
where the surface slopes are large, or for large angles of 
attack. 
For use in flutter work, the third term in the pressure 
expansion has usually been found to be negligible and is 


dropped. Also, only terms which are linear in the dis- 
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Fic. 2. Experimental and theoretical lift and center of pressure 


in steady flow at JJ = 6.86 


placements are retained in computing the oscillatory 


forces so that, as used, piston theory is nonlinear in 
thickness but linear in displacement. 
attack e‘Tects can be studied provided the third term 
in the pressure expansion is retained, subject to the 
limitation that w a < 1.0. 

It must be emphasized that piston theory is more 
than an approximation of linear theory. This method 
is a legitimate high supersonic theory which agrees 
with the second-order theory of Van Dyke and is even 
more valuable due to its simplicity. 

Landahl has presented a high Mach Number theory!! 
which may be described as an inverse Rayleigh-Jatizen 
method—that is, as an expansion in inverse powers of 
JJ. An iteration procedure is used, similar to Van 
Dyke’s, and the solution is obtained as an expansion 
in powers of 1/./ and thickness. An interesting result 
obtained by Landahl is that second-order piston theory, 
which had been based on one-dimensional ideas, can be 
obtained from the more general potential flow equation 
to second order and for large /. 

Tangent-Wedge Approximation—The tangent-wedge 
approximation, discussed by Lees in reference 12, has 
been used with success to predict pressures for steady 
hypersonic flow. A basic assumption is that the shock 
wave is very close to the surface but not coincident 
with it. The pressure at any point on the surface is 
taken as being the same as the pressure ona semi- 
infinite wedge at the same flight Mach Number whose 
half angle is the local inclination of the surface. The 
derivation, as carried out by Lees, uses the oblique 
shock relations with assumptions of large Mach Num- 
ber and small wedge angles. 

It would be of interest to examine some applications 
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of this method to unsteady flow. The local surface 
inclination becomes time dependent, making the pres 
sures a function of time. The resulting expression for 
the local pressure can be expanded as a power series 
in the downwash velocity and, to third order, is 


D— 


I / 4 p21") (w a) + ty + 1) t] (w/a)? + 


Oe i Sin ' 
-| “ed { ) 


ily fF 4)?/s 
where wa < 1.0. 

Eq. (6) is identical with the piston theory relation for 
the pressure at a point, Eq. (5), up to and including 
terms of order (w/a)*. The deviation between the ex 
pressions is only 10 per cent in the coefficient of the 
(w a)*® term. 
culations only terms to order (w a)? would seem to be 


As previously mentioned, in flutter cal 


ordinarily required and the tangent-wedge approach 
would be identical with piston theory. 

It is of interest that a similar approximation for three 
dimensional slender bodies is possible, referred to as 
the tangent-cone approximation by Lees in reference 
12. Anapplication to a cone-shaped body will be given 


subsequently. 


Newtonian Approximation 


The Newtonian approximation may be used to pre 
dict pressures on an airfoil or body of revolution. In 
this representation for hypersonic flow, the fluid par 
ticles ahead of the body are undisturbed. Upon 
striking the surface, these particles lose their momen- 
tum normal to the surface, then flow tangentially with 
no further change in velocity. This concept implies 
that the shock wave coincides with the body surface 
while the flow is compressed into an infinitesimal layer 
sometimes called the hypersonic boundary layer, be 
tween the surface and the shock. This does not account 
for the change in the effective shape of the body which 
can be caused by thick boundary layers and viscous 
effects. The pressure relation is derived in reference 


14 for a two-dimensional profile and is given by 


This pressure is determined by the local change in fluid 
momentum, given by the first term in Eq. (7), with a 
correction added to account for the centrifugal force re 
quired to make the fluid follow the surface contour, 
the second term of Eq. (7). 

The function f which appears in Eq. (7) describes 
the surface of the airfoil and varies with time as well 
as position. For a pitching and translating airfoil 
this function is 


f = + [h(t) + (x — cx) a(t)] + 


Y(x) 


) 


jupper surface 
\lower surface 
The instantaneous slope of this surface, f’, multiplied 
by free-stream velocity, corresponds to the downwash 
velocity w appearing in Eqs. (5) and (6). 

Any point of the surface in shadow, hidden from the 
free stream by a portion of the airfoil, is assumed to 
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region of near zero pressure. As reference 12 


points out, viscous effects would probably increase the 
est due to boundary-layer breakaway. 
Forces and moments determined from Eq. (7) will be 
independent of Mach Number. Reference 12 points 
ut that this is characteristic of the Newtonian approxi 
mation and is true in fact for slender bodies at J > 10 
ind at lower Mach The 
tangent-wedge and piston theory methods also become 


Numbers for blunt bodies. 
essentially independent of ./ for 1J > 1, as can be 
seen by rewriting Eq. (6) as 


C 2} (1/M) (w/V) + 


me ane ] ) 


{\ 
M[(y + 1)2/32] (w V)3} (9) 
ind noting that the first term vanishes as .1J becomes 
large 

An expression similar to Eq. (7) can be determined 
for bodies of revolution (reference 14+) which allows the 
Newtonian approximation to be applied to three- 
dimensional bodies. 


APPLICATION OF HIGH MACH NUMBER THEORIES 


Comparison With Experiment at Zero Frequency 


With several methods such as those just discussed 
available for use in flutter calculations, the question of 
which is best naturally arises. Some steady-state 
measurements of lift and moment have been reported 
ind can serve as a limited basis for comparison. 

In reference 15, McLellan, Bertram, and Moore have 
reported the lift and center of pressure as measured 
on some 5 per cent thick airfoil sections at .1/ 6.86. 
rheir results for a double-wedge airfoil are shown in 
Fig. 2 ¢ 
with predicted lift and center of pressure by the vari 


is a function of angle of attack and compared 


ous methods. 
The predicted lift coefficients are in good agreement 
with experiment by all methods up to angles of attack 
The predicted lift is identical by the 
Dyke, and nonlinear-pressure methods, 


of about 10 
linear, Van 
while Landahl’s technique is identical with second 
rder piston theory. The best agreement with experi 
ment actually results when third-order piston theory 
is used. 

rhe results of using a nonlinear thickness theory be- 
come more apparent when the center of pressure or the 
moment is examined. As Fig. linear two 
dimensional theory predicts that the center of pressure 
All the thickness 


2 shows, 


will be at the 50 per cent chord. 
theories show that the center of pressure will be shifted 
forward about 10 per cent of the chord from the linear 
theory prediction. The experiment shows much better 
agreement with the nonlinear theories, even to large 
angles of attack. This would indicate that these non- 
linear methods should give more accurate moments, 
even at angle of attack, and accurate moment predic- 
tion is well known to be of vital importance in aero- 
elastic calculations such as flutter. 

Van Dyke’ presents a limited amount of experimental 
data at \J = 1.42 and 1.61 on an oscillating biconvex 
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These data show 


airfoil which was obtained at the 
Laboratory, Teddington, England. 
a marked improvement, even at such low Mach Num 
bers, of the Van Dyke nonlinear theory over the linear 


theory. 
Comparison of Theories Through Flutter Calculations 


The various aerodynamic theories may be compared 


For this purpose, a 


by their prediction of flutter. 
simple two degree of freedom system similar to Fig. | 
was used with a 6 per cent biconvex airfoil. Flutter 
calculations were made over the Mach Number range 
1.3 to 10 for the elastic axis at 40 per cent chord, the 
center of gravity at 50 per cent chord, and a bending 
The 


sented in Fig. 3 by use of a parameter found convenient 


to-torsion frequency ratio 0.5. results are pre 


in comparing flutter results, (bw, @)Vu, which is 
called the stiffness-altitude parameter. 

The stiffness-altitude parameter depends only upon 
the physical properties of the wing —in particular, the 
torsional stiffness——and upon the atmosphere in which 
Its value increases as either altitude or 
stiffness increase. When _ plotted 
against Mach Number, as in Fig. 3, constant dynamic 
pressure curves will appear as radial lines through the 


The stable region, when this parameter is the 


it operates. 
as the ordinate 


origin. 
ordinate, will be above the flutter boundary. 

The most prominent feature of Fig. 3 is the diver 
gence between the flutter boundaries as predicted by the 
linear and the nonlinear theories at the higher Mach 


Numbers. For example, at Mach Number 6, the non 
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“a 
Fic. 4. Flutter boundary variation with frequency ratio for a 
biconvex airfoil, 1/1 = 5, r = 0.06, wu = 250, ra? = 0.25 


linear theories require about SO per cent more stiffness 
than linear theory in order to avoid flutter at a given 
altitude. 

For Mach Numbers greater than 2 or 3, flutter calcu- 
lations by Van Dyke’s method (including terms to the 
cube of frequency), Landahl’s method (using three 
terms), and piston theory give essentially the same 
The flutter boundary by the Landahl method 
In the lower supersonic 


result. 
falls between the other two. 
Mach Number range, the results of both piston theory 
and Landahl’s method deviate considerably from the 
Van Dyke boundary. However, the Landahl method 
has been shown to give better agreement with Van 
Dyke down to JJ ~ 1.5 if more terms are used in the 
expansion.'! The nonlinear-pressure results, using terms 
through the seventh power of frequency, fall very close 
to the Van Dyke result to JJ ~ 8, indicating that 
most of the nonlinearity due to thickness is contained 
in the method. 

An interesting point is the crossover of the flutter 
boundaries near JJ = 1.5 as predicted by linear theory 
Below this Mach Number, 
At higher Mach 


and the nonlinear theories. 
the effect of thickness is stabilizing. 
Numbers, thickness is increasingly destabilizing. 

In reference 16, for example, it is pointed out that 
the thickness effect varied drastically with frequency 
ratio when piston theory was used for the flutter calcu- 
lation. The other nonlinear aerodynamic theories 
predict this same variation as shown in Fig. 4 where 
stiffness-altitude parameter is plotted against fre- 
quency ratio at J = 5 for the same parameters used 
At this Mach Number, Van Dyke and piston 
The nonlinear 


in Fig. 3. 
theory give almost identical curves. 
pressure tnethod agrees essentially with them also. 
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The thickness effect is destabilizing at frequency ratios 
less than about 1.0 and stabilizing for higher frequency 
ratios. 

Flutter as calculated by the Van Dyke, Landa 
and piston-theory methods shows agreement for Jf > 
2.0 to 3.0. ’ 
latter two methods are not expected to give reliabk 
answers. In this range, 1.0 < JJ S 2.0 to 3.0, Vay 
Dyke's nonlinear thickness theory is probably the ap 
propriate procedure to use, provided cubic frequency 


\t lower supersonic Mach Numbers, thy 


terms are sufficient. The nonlinear-pressure method 
contains most of the nonlinearity but still does not give 
quite as good agreement with the more exact Van Dyk 
method for JJ = 2.0 to 3.0 as the much simpler theo 
ries. In the lower Mach Number region, it shows better 
agreement with Van Dyke. 

In the high Mach Number range, 1/ 2 2.0 to 3.0 
piston theory flutter calculations show very good agree 
with calculations made with the more exact 
They show better agreement than calculations 
In view of 


ment 
theory. 
made by the nonlinear-pressure method. 
its simplicity and in conjunction with its agreement 
with more exact methods, piston theory will now be 
used to show calculated flutter trends at high Mach 


Numbers. 


Flutter Trends at High Mach Number 


The nonlinear thickness theories just discussed have 
been shown to predict changes in flutter boundaries at 
high Mach Numbers when compared to linear theory 
The parameters used in Figs. 3 and 4 were deliberately 
chosen to emphasize the large variations which are pos 
sible. With a change of parameters, the thickness 
effects may be much less or the trends may be altered 
With the aid of piston theory, these changes due to 
thickness may be examined. At the high Mach Num- 
bers being considered, first-order or zero-thickness pis 
ton theory will predict the same flutter speed as linear 
theory and will be used for comparison. 

Frequency Ratio—Examples of the magnitude of the 
thickness effects for various parameters, as well as 
the change in trends with frequency ratio which may 
be expected, are shown in Fig. 5. A stiffness parameter 
is given at w = 150 for a zero thickness plate and for 3 
and 6 per cent biconvex airfoils. Fig. 5(a) is for the same 
elastic axis at 


) 


set of parameters used in Figs. 3 and 4, 
10 per cent chord and center of gravity at 50 per cent 
chord. At very low frequency ratios, zero thickness 
theory shows no flutter while airfoils with thickness 
have a finite flutter speed. However, above frequency 
ratios of about 1.2, thickness becomes stabilizing. 
Fig. 5(b) is for a 50 per cent elastic axis with a 60 
per cent center of gravity. In this case, the zero thick- 
ness wing has a finite flutter speed at low frequency 
ratios. The effect of thickness is large, changing from 
destabilizing at low frequency ratios to stabilizing for 
w,/@_ 2 1.1. Fig. 5(e) shows a case where thickness is 
destabilizing for all frequency ratios. 
here place the elastic axis at 40 per cent and center of 
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gravity at 60 per cent chord. Finally, Fig. 5(d) gives 
an instance with very little thickness effect at low fre- 
quency I itios. This is an ait elastic axis configuration, 
at 60 per cent chord, with the center of gravity at SO per 
cent chord. The crossover with frequency ratio is 
again apparent with a large stabilizing effect at ratios 
greater than }.0. 

Center of Gravity and Elastic Axis 


that thickness effects on flutter speed may have large 


Fig. 5 has shown 


variations which depend on the elastic axis and center 


IGH 


MACH NUMBERS ve 
variation with these parameters are shown in Figs. 6 
and 7 at V/ 
ratio of 0.5. 


5 for a bending-to-torsion frequency 
In both Figures, the stiffness-altitude 
parameter for a thick 
compared with that of a zero thickness flat plate. In 


3 per cent biconvex airfoil 1s 
Fig. 6, which is plotted against axis position, a thick 
ness of only 3 per cent is seen to require that the elastic 
axis be placed about 5 per cent farther forward to avoid 


Fig. ‘ 


flutter with a fixed amount of static unbalance. 
gives the flutter boundary variation with center ol 
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Fic. 6. Flutter boundary variation with elastic axis location 


for a biconvex airfoil, AJ = 5, wp = 250, re? = 0.25, wr/wa = 0.5 


a given axis location, thickness requires the center of 
That is, 
more mass balance will be required by thickness. This 


gravity be farther forward to avoid flutter. 


Figure also shows that, for a fixed elastic axis, there is a 
center of gravity location which gives a minimum flutter 
speed. The two Figures show that, generally, the 
effect of thickness is to increase the region of instability. 

Airfoil Shape—The previous examples illustrating 
effects which may result from thickness were for a bi- 
convex airfoil. The airfoil shape itself may influence 
the flutter boundaries, as shown in Fig. 8. Here a 
stiffness parameter is given for three 4+ per cent thick 
airfoil sections and a zero thickness plate. The 65A 
and double-wedge airfoils, when compared to the zero 
thickness plate, present curves similar to the biconvex 
airfoils. A large reduction in flutter speed is apparent 
at low frequency ratios with a crossover near w/a, = 
1.2. The shape effect is stabilizing at higher frequency 
ratios. The wedge presents a different picture—the 
result of its shape is destabilizing at all frequency ratios. 
However, the magnitude of this flutter speed reduction 
is small over the entire range. Generally, a blunt 
trailing edge has been found to decrease the size of the 
thickness effect on an airfoil. 


Possible Refinements 


At extremely high Mach Numbers, aerodynamic 
forces are altered by effects such as strong shock 
waves, viscosity, dissociation, ionization, etc. The pre- 
viously discussed methods do not account for such 
features in the flow. However, piston theory presents 
some interesting possibilities for use at these Mach 
Numbers, several of which will now be discussed. 

Local Flow—In the original work of Hayes, in which 
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he reasoned from the assumptions of slender bodies 
and large Mach Numbers, piston theory was based on 
free-stream conditions. In the physical picture as 
used by Lighthill, each point on the airfoil is regarded 
as a small piston moving in a one-dimensional channel] 
The pressure on this piston is determined by its own mo 
tion and by the speed of sound and density at its face 
which are free stream in the usual applications. This 
pressure is independent of the pressure at any other 
point. Thus, each tiny “‘piston”’ on an airfoil surface 
will feel only the local conditions and will not be in 
fluenced by free stream. Since Mach Number, density 
and velocity all vary locally, the pressure distribution 
could be altered and change the total force and moment 

Forces and moments as determined using the local 
flow conditions with piston theory have been used in a 
flutter calculation for a 4 per cent thick double-wedge 
airfoil. Local Mach Number, density, and velocity 
were calculated by the standard shock-wave relations, 
The results are plotted in Fig. 9 as the stiffness-alti- 
tude parameter versus flight Mach Number and com- 
pared with standard piston theory calculations. The 
local flow effect is stabilizing and is generally small. 

The local flow conditions used with piston theory 
also present possibilities for studying angle of attack 
effects. 

Blunt-Nosed Airfoils 
for the piston theory expansion in powers of w a to give 


Lighthill has pointed out that 


the pressure at a point, wa must be less than 1.0 
where w is the downwash velocity at the surface of the 
airfoil, given by Eq. (4), and a is the free-stream sound 
speed. If the theory is limited to small displacements, 
the downwash from Eq. (4+) will be approximately 


wm~ VY'(x) (10 
Therefore, w/a = MY'(x) < 1.0 (11 


is the requirement for piston theory validity. This 
shows that pressures will be inaccurate in regions where 
the airfoil slope is large. 

Airfoils with blunt leading edges are of great practical 
significance due to heating requirements. Over these 
leading edges, the slope of the contour is very large and 
the piston theory requirement of Eq. (11) 1s not met 
The Newtonian approximation, however, can be used 
to obtain pressures in this region of the airfoil, while 
piston theory can be used over the remainder of the 
surface. The effect of the detached shock wave should 
also be included for a more complete picture but has not 
been attempted in this discussion. 

This technique has been applied to a blunt-nosed 
airfoil whose surface is parabolic over the leading 20 
per cent and biconvex over the remainder of the chord. 
The airfoil contour is described by the equation 
V(x) = 0.84 reV xc, 0 <x < 0.2 

= re [0.117 + 1.64 (vc) — 1.76 (x/c)*], 12 


0265 x5 A 


Over the leading edge, the slope of this contour will be 
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Y'(x) 0.42 rV c/x 
The piston theory requirement of Eq. (11) will be satis- 
fied if 

Ww Y’(x) (C.42/V x) (Mr) Ve < 1.0 

Vx> 0.42 (Mr) Vc (13) 
The airfoil selected for the calculation was 3 per cent 
thick. Over the forward 5 per cent of the airfoil, New- 
tonian theory, as given by Eq. (7) without the centrif- 
Piston theory was 
The maximum 


ugal force correction, was applied. 
used on the remainder of the surface. 
slope occurring with the piston theory will be at x = 
(.05c, so that the piston theory aerodynamics should be 
valid for WS 17.7. 

The results of this calculation are shown in Fig. 
i0 and compared with a conventional piston theory 
calculation. For the frequency ratios below 1.0, re- 
vising the leading-edge pressures according to New 
tonian theory is seen to have a minor stabilizing in- 
fluence over a wide Mach Number range. When the 
frequency ratio is greater than 1.0, very little effect is 
seen at the lower Mach Numbers, but for JJ 2 15, a 
destabilizing influence is seen at w, w, 1.2, while a 
highly stabilizing effect appears at a frequency ratio 
of 1.5 

An interesting secondary result appears in Fig. 10. 
For the higher frequency ratios, a, w, 1.2 and 1.5, 
the curves have a maximum at JJ ~ 10-15. These 
maximums would be design poiats similar to the tran- 
sonic maximums determining the stiffness required to 
avoid flutter at a given altitude for this configuration. 


Flutter of a Cone 


There are certain instances when a cone may be in 
stalled on aircraft—for instance, as an angle of attack 
indicator. 


bilities of this configuration at high speeds. 


The question arises as to the flutter possi 


In deriving the air forces for such a study, a steady- 
state theory due to von Karman’? was used in a manner 
similar to the tangent-wedge or tangent-cone approxi- 
mation. A philosophy similar to that of piston theory 
was followed in the use of the instantaneous down 


wash. The pressure coefficient as given by von Kar 


man is 

c (2/M*)K.? log (2/K,) (14) 
where K w a and K, is required to be less than 1.0. 
Performing a Taylor expansion about A, = 1, insertion 


of the appropriate boundary conditions for a rotating 
and plunging cone, linearization with regard to the 
displacements, and integration over the cone result in 
the following expressions for the oscillating lift and 
moment on the cone: 

Ly Q, Lo = rA S.\Vk, 

i, TA/SMk?, Ls = (1A/6Mk) — 2x Lo, 


Meo = (2xA/6Mk) — 2xo Lo, 


M, = 0, 
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My = (wA/6Mk*) — 2x) L 
Wy = (wA/4Mk) — 2xo [Me + Ly + 2x LF: 
where 
A = (1/3) + M@ log 2 — (M6)? + (1/6) (1/6) 


This aerodynamic input has been used in a sample 
flutter calculation for a 15° conical shell. The require 
ment for A, < 1 limits the usefulness of these calcula 
tions to 1/6 < 1 or M & 7.65. 
where a stiffness-altitude parameter has 


The flutter re 


The results are shown 
on Fig. 11, 
been plotted against Mach Number. 
gion is below the curve. Results are presented for 
three axes of rotation, x» 1/3, 1/2, and 2/3, where 
the critical Mach Number varies from 5 at a 1/2 
to 3 at x% 1 3. Note that the axis placed at mid 
chord represents the lowest flutter speed. Axis loca 
tions ahead of or behind the midchord require less 
stiffness to avoid flutter or, conversely, could be flown 


at a lower altitude without flutter. 


CONCLUDING REMARKS 


Some of the theories for two-dimensional oscillatory 
air forces for use in flutter calculations at high Mach 
Numbers have been discussed. These include linear 
theory, Van Dyke’s second-order theory, piston theory 
Landahl’s method, tangent-wedge and tangent-cone 
approximations, Newtonian theory, and a method 
utilizing the nonlinear-pressure coefficient. The meth 
ods which include nonlinear thickness terms give much 
better agreement with experiment at .\/ G.S6 for 
static lift and center of pressure. The methods in 
cluding nonlinear thickness terms have also been shown 
to predict flutter speeds which may be very different 
from the prediction of linear theory. 

The second-order solution of Van Dyke is based on 
the general potential flow equations. This method 
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corresponds to linear theory as Busemann’s second- 
order solution corresponds to Ackeret’s linearized solu- 
tion for the steady case. As currently available, the 
method is limited in the lower supersonic region (J < 
2()) by having only terms to cube of frequency. 

For Mach Number greater than 2.0 or 2.5, piston 
theory may be used in flutter calculations. It has been 
shown to agree with calculations by the more exact 
method of Van Dyke at these higher Mach Numbers. 
Landahl’s method gives slightly better agreement than 
piston theory with Van Dyke when used in flutter calcu- 
lations. Agreement of piston theory with the tangent- 
wedge approximation shows that this simple method 
should also be valuable in the hypersonic region. 
Newtonian theory has possible applications at hyper- 
sonic speeds and around blunt airfoils. 

The nonlinear pressure method has been presented 
as an easy way to obtain an effect of thickness with a 
large number of frequency terms. A second possible 
use for this method is in conjunction with linear three- 
dimensional theory to obtain an effect of thickness to- 
gether with aspect ratio. 

The simplified flutter calculations presented have 
shown that thickness effects on flutter speed may or 
may not be important. Depending on the parameters 
of the system, such as frequency ratio, elastic axis 
location, center of gravity, etc., the thickness effect 
may be stabilizing, negligible, or destabilizing; the 
f the airfoil has been shown to influence these 


shape 
results 

Some possible refinements on piston theory for use 
at extremely high speeds were discussed. In general, 
the effect on calculated flutter speeds was small. In 
the absence of available experimental data, their value 
cannot be assessed. 

The possibility ot flutter of a cone at hypersonic 
speeds was shown by using the tangent-cone approxi- 
mation for aerodynamic forces. The critical location 
for the elastic axis was discussed together with the exist- 
ence of a critical Mach Number for flutter. 

The thickness effects discussed here have been shown 
to be capable of producing sizable changes in flutter 
speed. Studies of aspect ratio when combined with 


thickness are now needed. Information is also needed 
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on such nonlinear effects as b uundary layer, viscosity, 
shock waves, etc., since intelligent design for flutter 
demands an accurate knowledge of the aerodynamics 


involved. 
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Compressive Buckling of a Long Simplv 
Supported Plate on an Elastic Foundation 


PAUL SEIDE* 
The Ramo-Wooldridge Corporation 


SUMMARY 


\ solution is obtained for the problem of the compressive 
buckling of infinitely long, flat, simply supported plates resting 
on an elastic foundation. The effect of nonattachment of the 
plate and foundation is found to reduce drastically the buckling 


load of the plate as compared to a plate with attached foundation 


SYMBOLS 


A, B = constants appearing in deflection functions 

b = width of plate 

( = fraction of half-buckle wavelength representing length 
of region of plate in which the foundation is com 
pressed 

D = flexural stiffness of plate |D = Et8/12(1 — v? 

; = Young’s modulus of plate material 

F(x) = function representing variation in x-direction of the 
plate deflection 

k = buckling load coefficient (k = b?.V/7?D) 

k = foundation modulus 

N = critical compressive load per inch of plate width 

t = plate thickness 

Ww = plate deflection 

x = distance along plate longitudinal axis [see Fig. 1(b 

y = distance along width of plate 

= foundation modulus parameter (y = b*k/7*D 

v = (0*/Ox') + 2(0*/Ox?0y?) + (0*/Oy 

d = half-wavelength of buckle 

I = Poisson’s ratio 

Subscripts 
l refers to portion of plate which compresses foundation 
2 refers to portion of plate which pulls away from 


foundation 


INTRODUCTION 


A CONSIDERABLE LITERATURE exists concerning the 
behavior of plates on elastic foundations. In 
most or all of these investigations, however, the plate 
is assumed to be attached in some way to the founda- 
tion and is found to be considerably strengthened 
thereby. 

The problem investigated in the present paper is that 
of the effect of nonattachment of the foundation on the 
compressive buckling load of a long plate, simply sup- 
ported along its longitudinal edges. The analysis pre- 
sented here indicates, as might be expected, that the 
unattached foundation raises the plate buckling load 
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only a relatively small amount as compared to that of 
the plate with attached foundation, for which the 
buckling load is limited only by the inelastic properties 


of the material. 


ANALYSIS 


We consider an infinitely long flat plate, simply sup- 
ported at its long edges, and compressed in the longi 
tudinal direction. The plate rests on, but is not at 
tached to, a foundation of the Winkler type—that is, 
the deflection at any point is proportional only to the 
load applied at that point [see Fig. I(a)]. Since the 
plate is infinitely long, it can be expected to buckle in 
some repetitive pattern so that one part of the pattern 
compresses the foundation while the remaining portion 
separates from the foundation |Fig. 1(b) The de- 
flections of the former are governed by the differential 


equation 


DV uw, + N(07w, Ox)" — kw 0, 


and of the latter by 


+ N(07w2/Ox2”) 0, 
—(l-cA SxS (1—c)A 


DV 2*wWe 


The deflections of the two regions are connected by 


the boundary conditions 


Wi 7, , Wo} , ; : 0 
Ow; OWs 

Ox |2, . Oxe2 |, \ 

"20 O° Ws > 
a OXe? |, : 
O* 2); O* zen 
ONE;* |... » Ox? |, 1 . 


implying zero deflection at the boundaries of the two 
regions (by definition of the extent of the two regions) 
and continuity of slope, moment, and shear. In addi- 
tion, the deflections of each region are symmetrical 
about x, = 0 and x. = 0, respectively. 
We assume deflection functions of the form 

w, = F,(x1) cos (ry/b)| 1) 


io = F (x2) cos (1ry b)f 


Substitution of these functions into Eqs. (1) and (2) 
and use of the assumed symmetry conditions yields 


\ 





Fy(a 1, cosh ay(x, b) cos ax(x, by = j 
B, sinh a;(x,/b) sin (aex. 6b (5a 
i 1» cos 8)(x2/b) + Be cos Bo(x2/b 
where a (x/2 V2(V1 + 7) ¥ (k — 2)] 
) (r 2) (Vk + VR — 4) | ob 
1 bik/m'D, k b?>N/ x?D 





Substitution of Eqs. (5a) into the boundary conditions 





3) vields 
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cosh aye(A b)cos axc(A/b) + Si 3 
, 3 
eae, B, sinh ayclA b sin aoct(X b QO (Oa i 3 | ; ; 
h the SATITTTITITTTTATITVIVIVITVIVISTTTT IV ITTV ITV ITT 
—— 1. cos 3,(1 ( hb) + Fic. 1 Notation for problem discussed i) JLo; C 

By cos Bx(1 — c) (Ab Q (6b scone cen Noy ies ond ro : orn oe 7 
lila sinh a;c(A bh) cos avc(X b) — as cosh a,c(A b sin aoc(A b T 

’ SUp- B, |a; cosh ayc(Ab) sin acc(A, 6) + ae sinh ayc(X bd) cos aec(A/b 
longi AoB8, sin B)(1 — c) (AE) + BeBe sin Be (1 c) (A/? Gc 
t at- 
at is, 1; a — @ CC sh ay A b COS avc(A b = Paya: sinh ayclA b sin At A b a | 
- the ) B,[(a;° “= Gs sinh ayc(X b sin aec(A b) + 2 aya cosh a,c(A 5) cos avc(A h 
e the — 128;" cos By(i — « \/b) — BB." cos Ba(1 ( Nb 6d 
<le in : 
ttern l,iayla 5a sinh ay A b COS asc(X b) + Qa _ Day cosh ay \ b sin aoc(X b T 
rtion B,lay(a;? — 8ay*) cosh ayc(A 0b) sin avc(A/b) — az(ar® — 3a,") sinh ayc(A 6) cos avc(Xb 
. de- — 108,° sin 6)(1 — c) (Ab) — Bob.* sin B.(1 — c) (AVb Ox« 
ntial } 


Eqs. (Ga) and (6b) yield relationships between .1; and B,, and -1, and Bz, respectively The remaining three equa 
tions yield three relationships between B, and B, which must be identical. These, then, are the stability criteria 


(1) sinh 2ayc(A 6b) 2a,| + [sin 2aec(A/b)/2ae cosh 2a,c(A/b) + cos Zavc(rX b 
6, tan p; l " (Xr b ae De tan B.(1 _— 5 \/ b) De —_ D1 
_ = 3a") [sinh 2ayc(A b) /2a,| — (Q@ 3a lsin Pasc(A b) 2a = 
‘ 
) 8° tan B{1 — c) (A/b) — 2B. tan Be tan 82(1 c) (A/b 
4 by which are two simultaneous equations for k and c asa function of y and A 6 
Eqs. (7) can be arranged in a form that is better for computation by expressing them first as follows: 
A rN : ; lsinh 2ayc(A b 2a,| + sin Parc(r\ b 2a 
6, tan @ — ( — $8.tan B(1l — « (Bo- — Bi° Sa 
b b cosh Pa,c(A/b) + cos 2aoc(rA b 
} 
nr ) 
Bb,’ tan B — ( — — §5* tan Bo(1 — c) 
2) D b 
. a fs = 3a)" [sinh 2aic(A b 2a;| — (a;- — dae [sin 2asc(A/b 2a» 
(Go° — 1° : Sb 
cosh 2a,c(A/b) + cos 2aec(rA/b 
We next solve Eqs. (S) for tan 6,(1 — c) (A/6) and tan 6.(1 — c) (A 6) to obtain 
. rn (Bo — Qs” i. a 3a," |sinh 2a,c(A b 2a; | Tr (82° es - ja sin 2aoc(A b 2a: 
two tan 3,(1 — « - Ya 
ons) b 8,[cosh 2ayc(X/b) + cos 2acc(rA/b 
ddi- Dn (8,2 — a? + 3a,") [sinh 2a;c(A/b)/2a;| + (8)? + a? — 3a’ [sin 2acc(X/b) / 2as 
rical tan 6.(1 — « 9b 
b 8. [cosh 2a;c(A/b) + cos 2avc(dr/b 
Then tan 6.(1 — c) (A/b) = tan }(6./8,) tan! [tan 8,(1 — c) (A/4)]; 10a 
31° — ay? + 3a,”) [sinh 2ayc(A/b)/2a;] + (6:2 + ai? — 3a”) [sin 2acc(A 5) 2a 
}) oT — 
Be[cosh 2ayc(A/b) + cos 2acc(d/d) | 
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Variation of buckling load coefficient with foundation 
modulus parameter 
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Comparison of buckling loads of plates attached and 
not attached to the foundation 
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Fic. 2(c). Comparison of buckle half-wavelengths and region 
proportions for plates attached and not attached to the founda- 
tion. 


Ny 


Fic. 3. Limiting buckle configuration for plate resting on a 
rigid foundation 
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TABLE | 
k A/b 
0 $00 1.00 0.50 
] 1.33 0.94 0.44 
3 } 62 0 88 0.38 
5 t.7o 0.85 0.34 
10 } 92 O.85 ). 29 
20) 5.07 0.82 0. 26 
30 6.12 0.81 0.22 
+) 5.16 ().82 0.2] 
50 5.19 0.82(5 ). 20 
60 5.20 0.84 0.19 
70 5.21 0.83(5 0.18(5 
8.) 5.22 0.84 0.18 
OV) 5.23 0.83(5 175 
100 5.24 0.84 0.17 
5.31 0.84 0.095 


1,000 OD.¢ 


which is the desired stability criterion of one equation 
for k as a function of y and (A/d)c. In the solution of 
this equation for a given value of y, a value of c(\/b) 
was chosen and & was varied until Eq. (10b) was satis- 
fied. This process was repeated for other values of 
c(A/6) until a minimum value of k was obtained. This 
value of & is then the buckling load and, if the corre- 
sponding values of c(A/b) and (1 — c) (A 8) [obtained 
from Eq. (9a)| are denoted by x and y, respectively, 


' 


c=x/(x+y), A/b=xt+y 11) 


In the performance of the calculations, there was some 
question as to where the quadrant in which the angle 
represented by the inverse tangent in Eq. (10b) lay. 
A criterion was introduced to aid in determining the 
quadrant—namely, that the deflection functions repre- 
sented by Eqs. (5a) and (5b) must not change sign 
within the regions in which they are valid—a change of 
sign would mean that the region in which the plate was 
assumed to pull away from the foundation would con- 
tain a region in which the plate compressed the founda- 
tion, and, similarly, the region in which the plate was 
assumed to compress the foundation would contain a 
region in which the plate pulled away from the founda- 
tion. The numerical calculations indicate that, sub- 
ject to this criterion, the angle represented by the in- 
verse tangent lies in the first quadrant near the critical 
values of (A/d)c and k. 


NUMERICAL RESULTS AND DISCUSSION 


The calculated results of the investigation are pre- 
sented in Table 1. 

The values of \/b and ¢ are, by the nature of the com- 
putation procedure, less accurate than the values of &. 
These results are shown graphically in Figs. 2(a), 
2(b), and 2(c), where the results for a long plate at- 
tached to the foundation are also given. These are 


k “a et ee 
Nb=1/Y1+ 7 (12) 
( = 05 


It will be noted from Table 1 that the values of k 
seem to be approaching some asymptotic value. The 
limit of k as y becomes infinite can readily be deter- 


(Continued on page 394) 
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A Simple Method of Matric Structural 


Analvsis: 


Part I]I—Analvsis of Flexible 


Frames and Stiffened Cvlindrical Shells 


BERTRAM KLEIN* 


Convair Astronautics, A Division of General Dynamics Corporation 


SUMMARY 


The method presented in Part II (reference 2) is extended to 


cover the case of stiffened evlinders. Treatable loadings include 
producing peak shear, such as in flexible frame problems, 


Various structural irregu 


those 


ind those involving thermal stresses 
ceounted for, such as variable frame moment of 


larities may be a 


tapered longerons, variable skin gage, cutouts, noncircular 


inertia 
cylindrical shells, et \ flexible frame problem is solved for 
illustration purposes 
SYMBOLS 
cross-sectional area 
parameters 
frame spacing 
I = Young’s modulus 
F,5,f = functions of ¢ or factor 
= shear modulus 
effective centroidal moment of inertia of frame 
cross-section 
} = integral 
/ = developed length of outer line of curved element 
= frame bending moment 
= axial force or applied load 
= shear flow or shear flow coefficient 
= frame radial force 
= radius 
s = stiffness 
l = frame tangential force 
= skin gage 
= tangential and radial displacements 
J = frame shear 
= axial displacement 
= axial coordinate 
= coordinate 
, = thermal coefficient of expansion or function of ¢ 
= angle subtended by curved element 
= parameter 
A = determinant or difference 
= parameter 
€ = strain or error function 
= angular rotation 
¢ = angle 
K, = parameters 
= parameter 
= parameter 
= axial stress 
7 = temperature 
¥ = angular coordinate 


Superscripts and Subscripts 
. = effective 


= modified 


Presented at the Missile Structures Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958 


* Head, Structures Theory and Methods Group 


0, 1,2, = numerical indices 
( = centroidal 
H = horizontal 
? = numerical indices 
VW = relating to bending moment 
mm, rr, ete. = denoting various stiffnesses 
= outer 
R, = radial 
Ra = tangential 
V = vertical 
Ar = due to temperature difference 


INTRODUCTION 


IT Parts I AND II (see references | and 2), for the 
structures considered, it has been assumed that the 
shear flow is constant along each nontapered shear panel 

i.e., this element takes shear only. This assumption 
is simple and convenient but not always realistic. In 
flexible frame problems, for example, it is known that 
the shear flow can vary radically circumferentially. 
Therefore, in the present paper a method is developed 
in which the shear flow is allowed to vary along each 
curved element. The variation is assumed to be in 
the form of a power series in the angle subtended by 
the curved element. However, any other convenient 
series could be used, such as a Fourier series. 

The structures considered here contain, in general, 
curved segmental elements (which may be of nonuni 
form cross section), curved panel elements, and longeron 
elements (which may be tapered). The curved panel 
elements are assumed to carry axial as well as shear 
load as a consequence of the variable shear flow assump 


tion. The axial load carrying thickness is designated 











‘ 


Curved segmental element geometry 
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Fic. 2. Curved segmental element loading 
7; the shear carrying thickness is designated ¢. Both 
thicknesses may vary within the element. The radius 
of curvature, 7,, is assumed to be constant within each 
curved element but may vary from element to element. 
The dimension, 7,, however, may vary within each 
curved element. 

Loadings which may be considered include concen- 
trated loads either in the plane of the frames or in the 
axial direction, various distributed loadings in the 
axial direction, and thermal loadings of the frames and 
longerons. It is assumed, however, that there are no 
hoop stresses between frames and that the longerons 
cannot take bending or torsion. 


(1) CURVED SEGMENTAL ELEMENT 


(a) Shear Flow 


Refer to Fig. 1. Positive shear flow acting on the 
curved segmental element and the origin for the angular 
coordinate, y, are taken as shown. The shear flow is 


expressed as 
gq=gqgtalt/¢) + alv/¢)? + galW/¢)?... CG) 


The actual shear flow acting at certain points defined 
by xu, Yo, vs, ..., Wn may be expressed in terms of the 


coefficients, go, Gi, J, Jar. ss as 
lay.) [1 Gh/e) Wh/e)? (i/o)?  ) Fao] 
Vr 1 (poy) (pe ¢)? (Ye ¢)8 q1 
"Uy; = ] (v3 ¢) (Wz yg)? (Wz ¢)? ee q2 
q3 
Vy 1 (vn/¢) (vn/¢)? (dn/¢)3 




















(la) 


(b) Equilibrium Equations 


Refer to Fig. 2. The three equations of equilibrium 
of a curved segmental element are expressed in terms 


of the forces shown as 


vertical: (7; + Tisiddfi + (Ri + Ris) + 


(Fiigi + Fry3gs +...) 0) 
horizontal: (7; — Tiys)fi—' + (—Ri + Rin) - ° 
(240 + Fueog2 + Fusgs t+... dr 0 i 
moment: (7.7; — Fesed a) + (VW; -— Wyss) - 
[38qo + (8/3)de + (8/5)ds +... Ifo" QO 


The quantities, f and F are certain functions of ¢. 
These functions are defined in the Appendix and listed 


in Table 1. 


(c) Force-Displacement Equations—Constant Section 


The three force-displacement equations for a curved 
segmental element with constant cross-sectional prop- 


erties are 


—(M, + Mia) + (6; — 0541) Snm Ft 


(v; Vi41)Smr t+ (Ui — Uiai)Smt — 
(Fing: + Fags +... )re? 0 
—(R; — Riss) + (0; + Oj41)Sem + 
(Vi — Via) Sep Ht (tg + Ui) Sr — 3 
(Frogo + Frege + Frigs +... )r%o = 0 
—(7T; + Ty) + (0; — 0441) Sim + 
(v Visa) Ser + (Ui — Ui41) Su — 
(Figs + Frigs +... )ro = 0 


The functions, Fy, Fr, and Fr are defined in the 
Appendix. Certain numerical values are listed in Table 
1. The positive sign convention for “, v7, and 6 is shown 
in Fig. 2. The stiffnesses, S, are the same as those de- 
fined by Eq. (9) of reference 2. The equations are re- 
peated here in the Appendix, and numerical values of S 
are listed in Table 2. 


(d) Force-Displacement Equations—Variable Section 


When the cross-sectional properties of the curved 
segmental element vary, it is convenient to write the 
force-displacement equations in the following form 
(see the Appendix for the derivation) : 


—(60; — Oi41) + Fy M, + fr T; a Fr R — [(7, go + (d, gi + (Fo.)1g2 + } + (f4,)1 0 
— (Winn — e,,,9i41) sin B — (0; — 0441 cos B) + Fy.M, + f7.T; + Fr,Ri -— 
((Fe,)29o + (£o,)091 + (Fo.Jog2 +... |] + (Faz): 0 t 
Vi41 Sin B — (u; — Ui41 cos B) + (72,01 — rc,, 0141 cos B) — Fy,M, — F7,T; — fr,Ri + | 
[(Fo)390 + (4o,)3g1 + (Fq)ag2 +... | + Ca, 0 


The integrals appearing in Eq. (4) are 


B e3 
Fy, { (r./El)dy’; fr, = re tu, — | (7.2 cos W’/EI)dy’ 


*8 


B 
Lr, = f (y.? sin y’ EI)dy’ = 1 ee Fr = rody a, | (Ar,? sin 2y’ 2EI)dyp’ 


0 
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. 90, = =14191755592 16652767416 «0925618358 = « 7312980430 04321904666 «5936592365 3377308976 
urved 55  - 1e 428148006 =—-16795526008 = -« 1245561683 = «9755219384 =—--.5 745996387 «5768398938 «3234671602 
we 60 1. 732050807 1931341414 ©1630536876 14307973372 07675479626 = 5577529705 =» 35075145183 
65 26144506920 24059341625 020958620253 14780652174 10040603246 «5361690886 2890906825 
6705 26 414213562 26120120659 023607383507 9. 2006490236 14223611481 °5243567877 2791553677 
70 26 747477419 26178587769 = «2648731855 = 24 97081726 0452956894 «5118128411 2685770007 
7 Je 732050807 2 0288065035 03299288427 30702155605 20142777227 04843523475 «2455102513 
(3) 80 5671281819 24386677720  .4058490830 64123512817 32524625863 «4553856216 92195719348 
85 lle 43005250 2 473241717 04938180105 13.40926486 7671376727 04184225881 .1903741516 
the 
able TABLE 1(b 
own Functions of ¢ for« = dX = | = 1 
de- 
> Te- es FR Pro PRA fn "33 fa 3 
of S 2006092739612 2001160496297 0004835024589 0747883891 .04153640232 .00001266162877 .000004181730686 
02457645461 .004642571485 901934292589 «145079254 .0830798568  .0001013661295 .000C0337d375017 
«05486659699 .91044796L381 004352751509 02240894812 01245611212 0003425248613 0001141298969 
20975885281 201857951206 6007739683721 2984630857 1659719425 «0008132832611 0002709033577 
ved 123547454 -02351873614 00979664596 23355497442 «1866454240 6001159175993 «0003860478763 
= 01525777094 202904104729 «01209613309 03725576426 «2072882453 4001591925096 «0005300575310 
22198755926 004183725953 01742337517 4465005095 2484356047 002758198129 «0009179431588 
02995350876 °05697349779 02372286340 05196170852 «2895390641 004393781659. 001461445449 
3915957389 20744555134 «03099615862 59243505795 33042350143  —.006582653962 002188035311 
04961232122 20942891455 203924485938 26646605615 03711110715 00941153776 2003125987040 
06151634706 01164798762 204847038927 07362244816 4115759475 01297235005 0004504599243 
(4) 7427591553 01410323366 «05867405649 08070351640 4517893159 01735267844 005753371746 
3349395143 01679497320 20698566802 08770012645 «4917216656 02265628565 © 007504129350 
6039713136 01972328612 «0820184771 09460264894 53513421477 «02898 364546 -0095891630€ 
67e5 1.121817268 02127609897 008346635995 e9801546751 5510253918 032564685553 01076756936 
1 «207058899 02288793362 0951587788 1.014008969 5706184100 03644251555 201204238006 
12386914821 02628822318 «1092757128 1.080840584 6095164209  .04514653687 01489926240 
1579164523 02992288877 «1243658122 12146406233 6480002781 05521590457 «01819703597 
8 1.783619387 03378986217 =o 1404255446 12210583054 6860520048  .0667780457 02197484977 
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73 B 
tr, = | (r.3/ EI) [(1/2)(« + d) sin? wy’ + (1/2)(k — A) cos? p’|dy’; fu, = [ (r.2 cos w’/ EI)dy’ 


0 J 0 


Fr, = retu, — (r.3/EI) [(1/2)(« + A) cos? y’ + (1/2) (x — A) sin? py’ ]dy’; 


e 


*p 
Fp, = | (Ar.3 sin 2y’,/2EI)dy’ (4a 


0 


where « = 1 + (El/r2EA) + (EI/r2EA*); X 1 + (EI/r2EA) — (EI/r.*EA*) fa’ 
(the coordinate y’ is measured from the end of the curved segmental element as shown in Fig. 1) and 
bd ° 
(4,1 = re | (r2/ED) } ("CAF n/e"] — (CUFF) na/e" 3] +... (-1)" (F of dy’ = 
0 
8 n 
P [ (r2/ED) > ["C—1)! (8) nae" Jd’ 
J 0 1=0 
eo n : 
(70,2 = re | (r3/ED) dD [(—1)! "C,/¢"—'] UF) nt — f(Hn-A(EI/r2EA)] sin py’ + tb 
0 1=0 


(F) n41—-1 (EI/r2EA*) cos py’ dy’ 


i — S(5),-; (EI/r2EA)] cos py’ — 


(F) nyi—(EI/re2 EA*) sin y’} dy’ 


eg n 
(f,,)3 = r | (7.*/ EI) > [(—1)! "C,/¢"—']  ((F), 


0 1=0 


where 
(F)o = sin yp’ (F"), = (&/r) 1")! / oe + 2] — 2S), 22 1 
(F)) = 1 — cosy’ (F’)o = (1o/re) (W’) — (Fc 
(Fo = (Y’) — (Fe 
(F)3 = (w’)? — 2(5), f =(n—1), (n-—-1) 21 th’ 
s if y’) 
(F)s (W’)? — 6(F)o f =1, (n—i)=0 
TABLE 2 


Values of Stiffnesses for Various Values of ¢g; x = A = 1 





degrees §_(r,/EI) sp(F,/£1) $4 (r,/E1) s_(re/81) (82/21) 8, (2/ EI) s,, (72/81) 
5 68.7300002_ -1967.826821 22492.36348 -392.5372409 44864721194  -772837 848 8833577 e02 

10 544632758669 = -490.5644322 2782.129145 -97-0861092 550602686 47724207817 2706566968 

15 = 22..84346318 + -216.9980171 8092847623 —42.37507963 158 «1459509 -9237 «68338 34475 50372 

20 =—-:17-08888514 = -121..2505965 33361332760 -23.22877101 63282052384 -2839731112 7802 «097109 

2265 15216651025 - 94143640 23023506516 -18.06346397 4360905971 17412577391 42042559758 

25 =: 13.62609504 = —- 7609541306 §=—- 1649857755 -14.36959537 30281569671  -1119.881922 2401 594553 

30 =: 11230916661 = =52..86196551 915596100 - 9456024881 16055883668 = - 5142933791 891891490 
964698756 -38 «34823928 54276696149 ~-62663592153 9 051659585 = 262.2006185 374.4612907 
8439397907 = =28 692934961 3447665633 —4.736944869 5 0704858745 = 143.2701465 17067427119 
7413727282 -22.47291014 22647291014 -3.505876787 3503876787 = 8242461850 82 22461850 
6.674365919 -17-85581931 1498281140 -2 6589822258 20173118901 = 48.83322057 40.97393737 
5973797999 -14.44089867 10.11162611-1.917308298 16342583744 -29.60650726 20075069956 
5427289177 -11.34481898 08386094281 «410040323 08140871605 =-18.11171771 10.45680221 
42960739997 9. 82575094 4 .581822911=1 2019452845 4753693419  =-11405618606 5 0146250102 
4.751976385 = 898075255 3,0719949510-0.85757 3284 03552184850 - 855922674 5 543347803 
4557120067 3 222502203 2 0993663196 = 071394138435 02598534129 = =_-6 6566125002 2 0397153460 
4,.203678617 —6 0935025796 16858234561 - «4721448963 01265108436 = = 30746481932 1003866808 
50891013654 -5.880702507 1036926517 -02791423903 204922033488 - 14920695745 03386704816 
32611896183 -5.2008407277 0438178859 =01243269561 -01087719923 - «747006988 0653546430 


(4b 


(4b’) 


f 
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and, finally 


tude [ (-. | EaArydA £1) dy’ 
J 0 JA 
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(F,,)o = { r, (; | EaArydA E1) ~— (| EaArdA EA) | sin y’dy’ 1c 
0 A JA 4 

(Fx. 3 = | r. | -(". | EaArydA £1) a (| EaArdA FA) cos y'dy’ 
0 1 A 


The last three integrals represent the effects of ‘‘thermal 


loading’ on the curved segmental element. The cross- 
sectional coordinate, y, is measured from the centroidal 
axis of the cross section positive when outward away 
from the center of curvature. The quantities, £/, 
EA, and -A* signify that E may vary with y in the 


expressions 


E] Ey*dA; EA = | FdA: EA* = | EdA* 
JA JA 1 
(4d) 


All the above integrals can be evaluated numerically 
on an automatic digital machine asa subproblem. The 
results can be stored in the machine in their proper 


matric location. 


(2) CURVED PANEL ELEMENT 


Refer to Fig. 3 for the positive sign convention for the 
shear flow acting on and the geometry of a curved panel 
element. Notice that as a consequence of the assump- 
tion of no hoop stresses the shear flow is constant in the 


y direction. 


(a) Equilibrium Equation 

The equilibrium equation relating shear flow and 
axial stress acting on the curved edges of a curved 
panel is 


—(7,0, — i-.o;) + (a r,) (dq dy) = (5) 


(7;/i;)o; + (2a/iL) X 
[qi + 2(¥/¢e)ge + 3(v/¢)*o3 +...) =O (5’) 


with the aid of Eq. (1). 


(b) Interelement Continuity Equations 


Consider the juncture of two curved panel elements. 
At such a juncture, it is specified that the axial stress 
and its first derivative be continuous. Therefore, the 
following derivatives must be the same for each element 


at the juncture: 


(1 /2r,) (dq dy), (1 /ir,*) (d*q dy*) (6) 


Evaluation of these derivatives using Eq. (1) yields, at 


the juncture, 


G1 + 2g2 + 393 +...) — | 
Pi(Qi — 2g2 + 3q3 — ... Jinn = O 
303 + 64 <6 se 


p2(g2 — 393 + 6gs — 


where 
= (L re (7 Bias); Pe L L 41)’ i beat 6b 


The above conditions are valid whether there is a longe 
ron at the juncture or not. However, when there is no 
longeron at the juncture, the additional requirement of 
shear flow equality must be specified—i.e., 

OG is ok i eee — 


(91 -— + eS “—~ ce 43 0) (¢) 
Notice that in the above equations the index / increases 
in the same direction as the coordinate y. 


(c) Boundary Condition Equations 


In addition to the above interelement equations, 
certain boundary condition equations must be satisfied 
as follows: (a) at an axis of symmetry, 

q = 0 and d*q dy? = () Sa) 
and (b) at an axis of antisymmetry, 
dq dy 0 Sb 
These equations may be evaluated at given points by 
use of Eq. (1). 


(d) Axial Stress-Displacement Equation—No Longeron 


This equation is (see Appendix 








a* 


Curved panel element geometry and loading 
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Fic. 6. Matric equation for cantilever cylinder under symmetric radial load 
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i If 7is constant and no temperature gradient 1s present, (e) Shear Force-Displacement Equation 
the equation 15 The average shear angle (times twice the panel area 
of the curved panel is* 
































Oj w, — w;) + (a 2E) (o;, + o;) = O Qa) 
Uy 
u Wisr + Wii) la + [Visa +2; — (8; H+ Vj41) |Frorm, + [ai + Uinn — (4; +H for, + 
a F *Y ®y 
iteicn A, + 6541 — (0; + 6524) fy (Br iG F 1 thdy + g tiv ¢jdy 
JoL Jy . 
_| ta vy } 
\ | q | (1/t) (W ¢tdy +q lt) (by ¢)* dy 4 dx (10 
oa J 
\ . . - ~ . 
7 If fis constant, the right-hand side of Eq. (10) is 
w 
IT /tG) flan + (1/3)ao + (1/3)a, + (10: - . ; — 
wal. IG) |q (1/3)g2 + (1/9)a 7... 10a) rhe necessary 39 equations consist of 
The functi “ro, fe, % f; are defined 1 wv. ndix. ‘ 
oi } The functions Fo, fs, and /; are defined in the ~~ lix 3 X 3 = 9curved segmental element equilibrium 
a re ‘alues are tabulated in Table 1. These Doe 
Numerical values re tabulated in Ta le | eSe } axial equilibrium 
runcti . é > que > de- »s tie 
functions are proportion il to the juantities a le 6 interelement continuity 
1locogu - . : fere re 2 1e nme : 
jue rived in reference 5—1.€., } boundary conditions 
Ar Fror, 2L\a,| = for, 2L?\a, fy? 3 X 3 = 9 curved segmental element force-displace 
(10b) ment 
\ ; ; 
' ! axial stress-displacement 
(2) LONGERON ELEMENT 3 curved panel shear-displacement 
Lad 
— pea ae ' ' The matric equation defining the problem may b 
(a) Equilibrium Equation _— t punts = g id Oo y be 
| written as shown in Fig. 6. <A particular structure is 
L— Refer to Fig. 4. The equation is analyzed having the following values of elastic and 
F —(4d,o; — A,0;) tal[—(g +a + gt... )it geometric constants: 
P| = oe ' } ) . - : 2¢)° 
| q q1 T q2 3 ~ 4) OG (11 E 10°. G | x 10°. io ws ws 30°. 
= - - : 7’ (2La tG) X 10 3.90 
(b) Force-Displacement Equation . : 
7 a . . ° » , ; ) V4 ) 
= The equation is (see Appendix) Y 10, a = S34, us Lt 2a) X 10 02.8 
+ “a i, = 9 25, k \ ;, 
[1 — (x/a)/EAldx + -) - 
* lio [1 x/a)/EA |dx 5 a 2E) X 10! 0.417 
LJ - ((v/a)/EAl] dx 4 The multipliers of po vers of ten for uw, 6, and y are used 
= 70 because the unit values of displacement and stress are 
. ; aken as 10~* 0%, respectively. tP = 2 
bi a le ah) ae taken as ind , respectively. Also let / od 
x J0 According to reference 4, the two basic parameters 
rq wr . defining the present type of problem of a concentrated 
=| If .1 is constant and no temperature gradient is present, . ing ¢ - I : YI I ree 
) ae = /On\. ; loaded flexible frame are 
z the equation is the same as Eq. (9a) above. 
x A numerical example is given below to illustrate the 4 6r°i/ Jas: B G6(G/E) (3/t) (r/a)? 13 
= workings of the method applied to a particular problem. 
= Numerical Exam ple—Consider a uniform cantilever The values of the elastic and geometric constants listed 
= circular cylinder with a uniform edge ring loaded sym- above have been chosen so that the values of the two 
+— metrically as shown in Fig. 5. Because of symmetry, 
4 only one-half of the structure need be considered. The * 
oO 
= elements are taken as three curved segmental and three 4 —REF4N=C=D20 — a | 
curved panel elements. Each element subtends an 12+ © THIS PAPER r \? 1 + 
4 angle of 60°. 39 unknown quantities are found to . = a cee ee ieee | 4 - . 
a exist as follows: er ‘a T T Pp T 
10 curved segmental element forces: 7), 1,, 72, JJ, Ob) | T ] | T 
Ca ae wis 04!- | | | | 
8 Rs, Ts, Ms, Ra, Ts, Ma aN 
13 shear flow coefficients: (Go)i, (Qi), (Gadi, (Gs)1; bail ae | 
= — 
Go)e, (Gud, (G2)2, (Gs)23 (Gods, (Gr)3 (Gz)s» (Ga)3. (Gs)s 0 40 60 ~80_ 100 120 140/160 180 
3 0 f1)2 2)2, (a)e; 2 i af — 1. DEGREES 
t axial stresses (at the fixed end of the cylinder): -04+— __\P, DEGREES  )—_,. —“ } 
: } J}, F2, 63, 4 -06 
12 displacements: 7, @}, t%2, V2, 02, We, U3, U3, 63, Ws : : 
j ee, ee ee Fic. 8. Moment distribution in the loaded ring of a one field 


cantilever monocoque for 4 = 104, B = 200 


a's 
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392 
parameters .! and B are nearly equal to 10* and 200, in Figs. 7 and 8. The circled points represent the 
respectively. Curves showing the shear flow and values calculated from the results obtained by solving 
bending moment distributions are presented in reference 
} for the cantilever cylinder with such values of 4 and 


the matric equation defined in Fig. 6. These results 
shown below have been calculated in less than half 4 


B. These curves are repeated here as the solid curves minute on the IBM 704 electronic digital machin 


T; —0.0273SS9S1 (qo): 0.007049858S9 a; = —O0. 15945223 
M, 0.98982704 — (gi)1 0.0058217902 a» 0. 25492525 
T. —0.14856292 — (q2); ~0.0061917007 o; = 0.32102969 W 
Ms 0. 16682440 (s)1 —().0020638406 — oa, —(0).952435S89 | . 
R, 0.088292236 — (qo)e —().019953208 ra —(0.0015006301 X& 10 
iz 0.50179642 (qi)2 — 0. 028488696 W —().00006649158S6 & 10 
M —2.7(444229 (q2)2 —0.0013518929 — 2% 0.00023987483 & 10 
R 0.019316998 — (qs)2 0.0036772564 — v» 0.0020238958 X 10 : 
T; 1 .0529797 (qo)s —().046835175 A. —0.00020120508S * 10 
M; S. 6105814 (q1)s 0.024666250 ws 0.000097962769 X 10 . 
(G2) 0.027570091 1 —0.0084783675 & 10 
(g3)3 = —0.0016135264 2 : 0.0027844979 X& 10 
(qs) —0.00378SS3632 6 0.00016191069 X& 10 
w 0.00013386935 X& 10 
c's —(.0086417254 X 10 
w; = —0.00039716752 X 10 
W 
The error function obtained by substituting the above results back into the matrix is It 
£1 —0.Q000007 1 €14 0. 0000002 1 €2 — 0. QO0000006 
€) —0.00000025 €15 = 0. 00000023 €3 = —O0.000000006 } 
: e n 
€ = 0. 000000015 e156 = —O0.Q00000000F e259 = —O0.0000008S9 
€4 0.0000000036 67 = —O0.00000044 €5) 0 
€5 —().QO00000008S8 — exs = 0. 00000017 €31 0. QOO00004 
al 
€ = —(). 00000059 €19 0. 00000025 €29 —Q. 00000024 
€7 —(0.00000019 €2) 0. QOO0000 1 €5 0 
€s —0.00000015 €1 = 0. 00000052 €34 0 
€9 0. 00000035 €09 0. OO000007 €35 = 0 Fe 
€39 = —0.000000070 €23 0. 00000018 €35 = —O0.000000002 
é1 = —0.000000047 €24 = 0 637 = 0. QO000007 eG 
é2 = 0. QO0000048 65 = 0. 00000012 €33 = —O0.00000004 ta 
€15 —().000000027 €25 0. 00000059 €39 0 .00000014 Fj 
. ‘ . ~ Ca 
The worst error is associated with Eq. (29). However, since only eight significant figures have been carried in the ot 
computations, the results may be considered very accurate as far as the matrix inversion is concerned. Examination tic 
of Figs. 7 and 8 reveals close agreement for the shear flow and very close agreement for the bending moment when 
compared with the values obtained in reference 4. 
APPENDIX 
The functions, f and F, appearing in various equations in the text, are defined as follows: 
f, = tang; fo = —([(Fro)e=r-1 — 211; fB=B—-f; fs = (tand/o) —1; fe (tan @/¢)~!' — 1 
Fy, = 2f4; Fv; = 4+ 611 — 26°) fs; Fue = 2(1 + 2677/5) 
Fis = 2[1 — 8p? + 46-2 (1 — 66>2)f5| . 
Fro = 2+ 2! (r,/re) [811 + cos 8) — 2 sin B] — (« — A) sin Bf A 
Fro = }2(r,/r-) 86~! [cos @ — (sin ¢/¢)] + 4[sin @ — (6/3) cos d}f /A,, + 
} 16A¢~' [cos — (sin @/¢)] — 4(x — A) (1 — 267?) sin o| /A,, — 2(1 — 26-2) — 4(xk — AO! (1 + B-'sin B) A, wi 
10 
Frs = S(r,/r-) | [1 — 126-2 (1 — 26-*)] sin @ — [(¢/5) — 46-! (1 — 66~2)] cos of ‘A, — el 
32d [36-2 (1 — 4672) sin @ — o7! (1 — 12>) cos g]/A,, — 4(«k — A) { [1 — 129-7 (1 — 267*)] X pe 
sin @ — 249? (1 + B-' sin B)} /A,’ — 2[1 — 126-2 (1 — 267° r 
} te 


Fm = (Fri + (Fre; Frs = (Fr3)t + (F73)2 





nt the 
solving 
results 
half a 
fr 
1¢ 
} 
) 
| 
) 
n the 
ition 
vhen | 
A, 
x 
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(Fri)i = S(ro/r-) (cos @ + [(g?/3) — 1] (sin ¢/¢)5 A,’ 
(Fr) = 8[— (x — 1) (sin Q >) + (1/2) &« - d) cos ¢] Ay + (26 


(Frs)1 = 4(r,/r-) | [(2/5)@? — 611 — 26-2)] (sin @/¢) + 211 — 66>?) cos of / A, 


F, 1}2(1 — 6@~*) (sin @/¢) — 6x[2¢~-* cos ¢ + (1 — 46>?) (sin @ g)] + 
(x — XA) (1 — 66-2) cos of /Ay + 66-' (1 — 26-? 
Fin = fs(Fr)i + (8/3); Fuss (r-/ro) Ufs( Fr3)2 — 267! Frys) + fs(Frs)1 + (8/5) (A-1 
where « = 1+ (/r.) (A + A®*); A = 14 (1/2) (A — A® (A-la 
and A, = «8 —dXsingp; A, A, sec @ (A-1b) 
A, = «8B + Asin B — (4/8) (1 — cosB); Ay = A, sec gd (A-Ic 


Numerical values of the various functions are listed in Table 1 for several values of », given in degrees. The func- 


tions involving x, A and (7, 7.) as parameters are presented only for the limiting case of « r (r,/¥7 g 
Stiffnesses 


The values of the stiffnesses appearing in Eq. (3) of the text are 


ae Sm: Sn EI rfaSalte?/ED) . 0 0 rfs hi FT| | l/¢) O O 
:, a ae = — 0 Sel hse / eel 0 —r.fy 1 fil t 0 0 O (A-2) 
 * Fs ‘ 0 0 Silt3/ED) | irs —h 1 ; 0 OO 
where S,-(r-°/EI) = 2(1 + cos B)/A,; Sxu(r-2/EL) = 2(1 + cos B)/A, A-2a) 
Itis to be noted that Sin = See and Sy = —fodrm A-2b) 
Numerical values are listed in Table 2 for various values of ¢, for the limiting case of x r 1. To obtain nu 
merical values for other values of x and \, multiply S,,,, S,,, and S,, by 
(8 — sin B)/(KkB — Asin B A-2¢ 
and the remaining stiffnesses by 
[8 + sin B — (4/8) (1 — cos B)]/ [kB + A sin B — (4/8) (1 — cos B) (A-2d 


Force-Displacement Equations for a General Curved Segmental Element 

rhe equations are derived by the use of the principle of virtual displacements.’ For each equation a system of self 
equilibrating dummy loads is required. The systems chosen consist of a unit positive moment, radial force, and 
tangential force applied in turn at the end of the element designated 7 (from which the coordinate ¥’ is measured in 
with the necessary equilibrating set of forces (and moment) applied at the other end of the element in each 
Internal work is comprised of bending, 


Mathematically, the three condi 


Fig. | 
case. The work done by the dummy loads in each system must vanish. 
stretching, shearing, and temperature distortion work of bending and stretching. 


tions become 


*3 %5 . 
0; — Or41 (r.W EI)dy’ + } r (| akArydA E1) dy’ (A-3a 
/70 70 JA 
— 041 COS B+ (Mist — Pe,,,0i41) Sin B = (CVUr,? sin ~’/ EI)dyv’ + (r-V cos py’ EA*)dy’ 
y /70 J 0 
} r.P sin py’) EA)dv’ + (44,)2 (A-3b 
~/7 
Ui — Ui41 COS B — BOiyi(%e, — Pe,,, COS B) — Viz41 Sin B = [r°\J(1 — cos wy’) /EI\dy’ + 
70 
bd] °3 
(r.V sin ~’/ EA*)dy’ + (r.P cos ~’/ EA)dy’ + r.(Fy,)1 + (44, (A-3e 
v7 JV 
Where JJ, I’, and P are the bending moment, shear 
force, and axial force acting in any cross section of the , ; ; 
1 , , ; : : 3 Force-Displacement Equation for a Longeron and Axial 
element, and the integrals are taken from point 7 to . , . . 
: Stress-Displacement Equation 


pont: + 1. It is convenient to add —r,, times Eq. 
A-3a) to Eq. (A-3c). In this manner Eqs. (4) of the Eqs. (9) and (12) of the text are derived directly from 
text are derived. the following basic relations: 
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Compressive Buckling of a Long Simply Supported Plate | 
Continued from pa a SA 
mined by noting that, when the foundation is rigid, the The value of ¢ is, of course, zero since the foundatiot 
required buckle pattern will be as shown in Fig. 5 does not deflect. 
The compressive buckling load of this configuration It can thus be seen that the lack of attachment of the 
is that of a flat plate with clamped loaded edges and plate to the foundation results in a drastic reduction in 
simply supported unloaded edges and is available in the buckling stress coefficient, from a value limite 
reference 1. The lowest value of k compatible with the only by the inelastic properties of the plate material 
criterion of no change in sign of the deflection function is to a maximum possible value of 16.5, or only a 55 pet 
ae a cent increase above the plate without a foundation 
k 16.5 la 
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and the corresponding value of Gis : ; 
Timoshenko, S., Theory of Elastic Stability, p. 864; MeGt1 
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Water-Based Helicopter Symposium 
Stevens Institute of Technology, Hoboken, N.J., June 27, 1958 
The first Water-Based Helicopter Symposium to be held at Stevens Institute of Technology is being sponsored 
a 


by the Bureau of Aeronautics, Department of the Navy, and Stevens 


Several papers dealing with the operational concepts and hydrodynamic requirements of water-based heli 


copters will be presented by representatives of the services, Government, industry, and university 
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On a Subsonic Compressibility Correction for 
Propellers and Rotors 


Svetopolk Pivko 
Aerotechnical Institute, Beograd-Zarkovo, Yugoslavia 
January 24, 1958 


fs EFFECT of compressibility for propellers can be approxi 
mately predicted by applying the blade-element theory, and a 
obtained in 


general subsonic compressibility correction can be 


the form of a ratio of the thrust and torque coefficient values in 
] 


subsonic flow to the corresponding values ia incompressible flow 


The effect of compressibility for each section of the propeller 
blade can be approached by the linearized Prandt!-Glauert cor 


As the 


ing inflow velocity distribution has a negligible effect ia predict 


rection for two-dimensional subsonic flow ictual vary 


ing the ratio of the propeller characteristics, the constant inflow 





can be assumed as given by the Theodorsen theory! and the mo 
mentum theory for compressible fluid.’ 
The thrust and torque of a blade element at radius ~ can be 
expressed as follows 
dT = (1/2)Bp, WC (cos @ sil @ tan e)c d l 
dQ = (1/2)Bp, 2 In S tan e¢ ( > 
where 
i | 2 2 = R*w*(r 
@ = tan! j [I 2 ; = (A/x) [1 2 \ 
Yr = R i = J 1 = V,/R \ A} 1 4 
In these expressions, B is the number of blades, p ir density 
it propeller, Rk propeller radius, c = blade-section chord, tan e 
= Cp/( ratio of section drag and lift coefficients, r 
sultant local velocity relative to the propeller blade element, 
issuing constant inflow and neglecting the slipstre rotation 
I) and w propeller axial advance and angular velocity, respe¢ 
tively, and x increase of the axial velocity in the Ipstream 
° 
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Substituting these expressions into Eqs. (1) and (2), and setting 


Cr (Cr l (W/a,)? (3 


where a, is the velocity of sound at the propeller, the values 7 and 


Q can be found by integrating Eqs. (1) and (2), assuming a con 


stant w, constant chord c, and constant lift coefficient for incom 


pressible flow (C, along the entire radius. For the correspond 


ing incompressible case, the integration of Eqs. (1) and (2) for 


the same propeller gives 7 and Q Thus, one obtains the 
ratios 
T/Tis Cr/(( = 

(p/p 1 B tan « A By tan « } 
O/O = Cy/(Ce = Cp/(Cy - 

(pi / p c D tan e ¢ Dy tan e 5 


», D, and Dy are functions of 


where 4, 4 


ind Wis Wie+ WV 


If one is justified in neglecting the term tan e = Cp/C 


ind (5) reduce to 


It is interesting to note that in this case the propell 


not dependent upon Mach Number, 


and velocity of sound a, at the 
wy = V, 


The first-order effects of density p 
propeller? may be easily obtained as functions of 


thus, 


10) 44 ; ' ’ ' ' r 1 
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Fic. 2 
si Iwi ; 
| 
a 
pPr/ f l O.50My*, l 0.1] 
where is approximately given by 
=w/V 7/2 1/2) [T/(p0/2)Vo?#R } 
rhe more exact magnitude of can be evaluated from t dat ) 
measured by Theodorsen! using an electrical-analogy 
Fig. 1 shows a graph for a four-bladed single-rotating pr 
from which the values of can be read directly witl 


laborious calculations. Since the effect of is very sm 


if One is interested in only the first-order compressibility correc 


tions, the substitutions 


WV Wy = 1] V/ = R 

and Mr My = (\? / 
where dy is the velocity of sound of undisturbed flow, may b 
introduced into Eq. (6 Fig. 2 shows a plot of Eq. (6) wit 
these substitutions, with p; = po and for constant values of / ; 
The dotted lines shown for comparison are those given by tl 
simple Prandtl-Glauert rule* 

C7 (C7 1/(1 VW Pi é 
using the helical Mach Number occurring at the three-quarter 
radius 

Ms, = , M,? + 3/4) Mp] *} 
From Fig. 2 it is evident that this simple compressibility corre igt 


tion provides a satisfactory approximation as long as the Mach 


Number remains less than 0.6 

Since the blade-section compressibility correction as given by 
Eq. (3) is restricted by the limitations of the linearized subsoni 
flow theory, it may be assumed that the corrections evaluated by 


Eqs. (4), (5), and (6) are limited by the critical Mach Number o! 


OS or 


t 


less. A comparison with available experimental data 


proves that the derived theoretical results are valid as long as 
the Mach Number of 0.8 is not exceeded 


results prove also that the compressibility has a negligible effect 


These experimental 


on the propeller efficiency [see Eq. (7 
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Analysis of a Transpiration-Cooled 
Hemisphere-Cylinder® 
k € na cal Laboratorie l e O M esota 
t] nalvsis of Stine and Wanlass! is pre 





EXTENSION of the an 
A sented, which predicts the effect of air injection on the heat 
sfer distribution over isothermal blunt-nosed bodies of revolu 
O dvance knowledge of the surface pressure distribu 
red. Results are presented for a transpiration-cooled 
im Mach Number of 3.6 


Mangler transformation and the Stewartson-I]lingwort! 


cvlinder at a tree-stre 


sformation were used to relate the problem in axisymmetri 


flow to two-dimensional compressible flow super 


nsional 





n to an approximately equivalent two-di1 





flow (* superscript) where the thermal solutions 
9 


wedge flow cases presented in reference 2 


pplied. The effects of the major assumptions—constant 


N ver, small temperature differences, and local ippli 
of the Falkner-Skan solutions—-have been discussed in 
ce 3. No additional assumptions were required to extend 
lysis to the case of air injection 


heat balanes ver an element of the surface, the 





ipplying these equations to a hemisphere-cylinder the local 


Mach Number and Euler Number (7) may be evaluated from 
idified Newtonian law, C,/C; = cos*6, if the flow just 
utside the boundary layer is taken to be isentropic. Once Eqs 
ind (5) are solved, the solution of Eq. (3) is found by a trial 
ud-error procedure using the tabulated solutions of reference 2 
rhe results presented in Fig. 1 display the distribution of the local 
it-transfer parameter along the isothermal surface. The uni 
decrease in local heat-transfer coefficients with air injection 


, Varies with 





ipparent. The value of the blowing parameter, / 
gular location and approaches asymptotically, as does the heat 


ransfer parameter, the corresponding flat-plate values 
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Note on Mechanical Behavior After Creep 


George Gerard 

Assistant Director of Research, College of Engine 
New York University 

February 21, 1958 


INTRODUCTION 


| es AIRCRAFT or missiles operating in a thermal environm«e 
it is of interest to consider the mechanical properties 
materials at elevated temperatures after creep has occurre 
For example, in the superposition of a sudden maneuver up 
steady-state loading, such data would be of interest 


COMMERCIALLY PURE ALUMINUM 


iracterist 


To obtain some information on the short-time ch 
of materials after prior creep, the following tests were conduct 
on 3003-0 aluminum (commercially pure The tests utiliz 
round tubes at 650°F 
conducted in three parts, as follows 
a) To establish the basic short-time stress-strain charact 
plied at a strain rate 


load was i} 


approximately 0.002 min after soaking at temperature for 


istics of 3003-0 at 650°F 


nt, 


of 
d 


on 


ics 
ed 
ed 


under compression or torsion and were 


er- 
of 
l 
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Fic. 2. Effect of prior shear creep on the short-time elevated 


temperature shear properties of 3003-0 aluminum after exposure 


at 650°F 


hour \t a strain of approximately 0.03 the specimen was un 
loaded and then reloaded to a strain of approximately 0.06 

b) In the first prior creep strain experiment, a specimen iden 
tical to that used in (a) was loaded to a stress level in the neigh 
borhood of the proportional limit after soaking at 650°F. for 1 
hour. The specimen was then allowed to creep under this con 
stant stress until a strain of approximately 0.05 had been reached 
The specimen was then unloaded and reloaded to a strain of ap 
proximately 0.06 in order to obtain the short-time stress-strain 
characteristics after prior strain 

c) The second prior strain experiment was identical to (b) ex 
cept that the stress level under creep conditions was roughly 
half way between the proportional limit and the 0.002 offset 
vield stress 

The results of the experiments on 8003-0 aluminum at 650°F 
are shown in Fig. 1 for compression loading and in Fig. 2 for 


shear loading. An examination of these Figures reveals the 


highly significant fact that for commercially pure aluminum at 
650°F. under exposure times up to 80 hours, the path of prio) 
creep straining has no significant influence upon the reloading path 
This is vividly demonstrated in Figs. 1 and 2, where all reload 
ing paths at approximately 0.03 strain, regardless of the path 
of prior straining, follow an elastic path until the original short 
time stress-strain curve is reached and then follow essentially the 


same plastic path. 


\LUMINUM ALLO\ 


Phe above results apply to a simple material which is highly 
suitable for fundamental creep studies. In structural applica 
tions, it is pertinent to consider the behavior of the more com 
plex alloys under prior strain conditions similar to those con 


sidered above Here the effects of aging and deterioration de 
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pend significantly upon exposure time and it is possible that t 
behavior after prior straining may differ from that observed 
the 8003-0 aluminum 
For this purpose, test data of Beck! on the tensile propert 
, 
i 


f Alclad 2024-T3 aluminum alloy sheet after prior ter 


at 5OO°F. are significant The results shown in Fig for 
hours exposure indicate a possible strengthening of tl 
due to aging. Those shown in Fig. 4 are definitely in tiv 
the deterioration of this material after 90 hours exposur 

It is significant to note that, in all cases shown in Fig 
$+, the short-time reloading path following prior cree] 
elastic for an appreciable stress range beyond that correspor 
to the creep stress level However, the deterioration in pr 


ties shown in Fig. 4 indicates that under sufficien p 





times, the short-time reloading path may not be elastic be 
the prior creep stress level 

To conclude, it appears that for 3003-0 aluminum t S 
time reloading path is independent of the prior strain 
the conditions of the experiments. For an alloy such (24 


under the test conditions of reference 1, the more restrictive 
clusion, that the short-time reloading path after prior str 
is elastic considerably bevond the creep stress le ( 
made However, there are indications that for an 
sufficient exposure times the short-time reloading pat 


be elastic bevond the prior creep stress level 
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Fic. 3. Effect of prior tensile creep on the short-time elevat 
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temperature tensile properties of Alela 
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alloy sheet after 20 hours’ exposure at 500°F 
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Fic. 4. Effect of prior tensile creep on the 


vated-temperature tensile properties of Alelad 2024-T3 alumunu 
alloy sheet after 90 hours’ exposure at 500°F 
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© that Note on the Lift of Slender Nose Shapes 
served According to Newtonian Theory 
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is independent of 
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coeihicient 








Fic. 1. Reduction of d wit plitude of init 
he problem is considered as a small 
) rturbat n the solution at zero incidence Phe following 
; B .. f pansion is assumed | | | 
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dad + ITC 1 l+a 0 (2 
where A ind IL are paramect« depending up thie ad.1 ter 
| ils of the calculation are presented in reference 2 properties, and geomet! # the colums 
While this last equation cannot be readily integrated for at 
REFERENCES trarv values of exponents, certai bservatior re pertine 
D s Jourr th Phe integrati f Eq. (2 t relati betw the ti 
\ 
: 5 7 ( e amplitude If defines the time corresponding 
x4 : Z creep ing as the value of ¢ whe + then for finit 
RA ( n. Pa | 70. I iOS 
I ee , , 
creep buckling times t ilues of the p its must be limi 
| 
. 
; 2. , , 
— 9 . “ 
. . : \ closer examination of the case » 2 and» 3 will reve 
On the Reduction of Maximum Loads in a 
. . . in additional restriction n this instance, the creep buckling 
Nonlinear Viscoelastic Columns* 
time will be 
p Harry H. Hiltor /2)t., = [(1 — A)/8] In [(4 4 
‘ Aeronautical Engineering. University of r $A /37) tan 2/a ! 
| and Consultant, Hughes Aircraft Comper 
é t where dp is the initial amplitude of the deflection of an equivalent 
reo y 94 258 elastic column obeving the stress-strain relation 
™ - : F e = (o/s o 5 
T E PI LEM of creep buckling of linear and nonlinear visc 
008 ¢ tic columns has been treated extensively in the litera ind having an initial curvature 
I described in th paper was performed under the yONSOF é = ywsine ) 
) t Deve pment La ratorie Huhe Aircraft ¢ npany 
iteful to them for perm nto] 1 this articl Following the same p 


rocedure as was used in obtaining Eq. (2 
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one finds that for the elastic column 
(A/no) [(1 — ao)"0 — (1 + ao)%] = yo — a (7 
For the present case, % = 2, Eq. (7) yields 
ao = yo/(1 — 2A), A = 2k(L/xn)? (P/A)’ (8) 


It is apparent that the maximum or critical value of the elastic 


load is achieved when 
Ac... @ 1/2 (ny = 2 (9 


Since this elastic case forms the initial conditions for the vis- 
coelastic problem, A in Eq. (4) must be less than 1/2, otherwise 
instantaneous elastic buckling would take place. A meaningful 


physical ¢,, must obviously be positive which implies that 


A noz 3 In{[(4 + ao?)/ao?]} ) l (10) 


(1 + }32 tan! (2/a 


Fig. 1 has been computed using Eqs. (8) to (10) and illustrates 
the reduction in maximum load for a viscoelastic column of this 


type 
A similar analysis for the cases nm = 3 and mp = O and 1, respec- 
tively, reveals no reduction in maximum viscoelastic loads. 


Hoff’ has shown that the elastic term in the stress-strain relation 
(1) seems of little importance in creep analysis and can be neg- 
lected (ky) = 0) after a moderately large amount of creep has 
taken place. However, it would appear that initial elasticity 
strongly affects stability considerations—i.e., the reduction of the 
maximum load. 

The present analysis further indicates that the mere use of a 
nonlinear viscoelastic stress-strain law, such as Eq. (1), gives no 
a priori guarantee that a finite creep buckling time can be 
achieved. Unless the inequalities (3) are satisfied, this particular 
stress-strain law will not be of any more value than linear visco- 
elasticity. It is difficult to say, at this point, whether or not 
other more nonlinear viscoelastic laws possess similar restrictions, 
because, in general, closed-form solutions such as Eq. (4) cannot 
be obtained. However, it must be remembered that when use 
is made of Eq. (1), the constants, k, ko, m, and mp) are determined 
by fitting this equation to experimental data. Such approxima- 
tions might introduce restrictions into the problem because of the 
aforementioned properties of the stress-strain law rather than 


because of inherent characteristics of the physical problem 
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On Heat Transfer in Slip Flow 


Stephen H. Maslen 
Lewis Flight Propulsion Laboratory 
February 24, 1958 


NACA, Cleveland, Ohio 


A NUMBER OF AUTHORS have considered the effect of slip on the 


heat transfer and skin friction in a laminar boundary layer 


* The reader is referred to reference 1 for a comprehensive bibliography 
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over a flat plate. Reference 1 considers this by a pertur 
on the usual laminar boundary-layer analysis while some other | 
studies—e.g., reference 2—have proceeded through consideratio; 
of the Rayleigh problem of the impulsive motion of an infinj 
plate. 

In all these cases, there appears to be an error in the calcula. 
tion of the heat transfer to the wall When there is slip, th 
transfer has two parts. One is the usual contribution due 


while there is a second part due to ‘‘sli 





conduction (kT, 


friction” (uuu, This term, heretofore omitted, represents t} 
work done directly on the wall by the shear This latter ter 
is nonzero only if there is some slip at the wall. Then the 


transfer, g, into the wall is 
= kT, T BUM, 


In reference 1, a solution is given for the slip flow over a flat 


plate. The solution can be found in closed form for almost 
condition of the gas provided only that the corresponding « 


tinuum solution isknown. Fora perfect gas having constant spe- 





cific heats but whose viscosity and conductivity are 
functions of temperature, the velocity and temperature profiles 


are given by 


t/u = Moi) — € Vil 7 i n dn } 

{ L J) 

= ‘ a yo 
ar ee = 7)(1 €z\ Ton| A + f(n)dr — F(» 

iia ie cae ea 
where 
n= (W/V2) V Baobka/?P 
3; = Uy 0) / pol T du dT F(n bu 


— 1) 1M2Pr uo + (RoT on/p 


Uo) | ¢ dn 
{ 


and where, for supersonic flow 


Modi V 7 
\ Re 6* x 
is the continuum displacement thickness, 


by a = V1r2 


have their usual definitions i 


In these expressions 6* 
and a, is related to the reflection coefficient, d, 
2 — a)/a M, y, Pr, and Re 
terms of free-stream conditions and distance along the plate 
In these equations, all the required boundary conditions have 
been satisfied except that on the wall temperature or heat trans- 
fer. Thus, the coefficient B of Eq. (3) is undetermined. Re- 


garding conditions far from the plate, the displacement eff 


of the boundary layer on the external stream has been allowed 
for and the profiles given go exponentially to the required limit 
for all values of B 

wall values of most interest cat 


Number 


From these expressions, the 
The 


using Eq. (1), 


found temperature, shear, and Nusselt 


1 1 = To + ec) (moB/k t+ 4iUV 4(M l x 
] ‘ 
[u CV To/b] Toni 9 =0 
C;V Rez = [(pu V Re 1/2)¢ = 
Zul tex|2V R 0 | ; 
Nur/V Re, = x(kT, + pun, )yoo/k(T, — T,,)V Re: = 


1 (ft. ¥. koTon + €.(0 


T, is the continuum value of the adiabatic wall temperatur 


is related to the accommodation coefficient, c, by 
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erturbatior f : , , , 
-— These relations follow from Eqs. (2) and (3 However, Eqs 30 


ee 5 wa be shown to be much more generally valid, requiring nei- -—--— METHOD OF CHARACTERISTICS 
msideratior 3 ‘ : 
‘ we 1] ther equilibrium nor any particular assumption regarding the APPROXIMATE_ METHOD 
= state equation or the transport properties. It has, however, 25 4 EXPERIMENT (REF. 2) 

; p r A , 


been assumed that the external stream is a perfect gas and that 


some other f 








he cal ula ‘ . _— 
— 1 forces and diffusion and radiation effects do not occur 


;M=200 6=10° 





is slip, this ga ; : 
; ; : Several results are at once apparent from Eqs. (4) and (5 
10On due : 


to “slid: If the thermal boundary condition specifies the wall temper- 
0 Siding . ’ . . 
a iture, then B is determined and the solution is completely de- 
resents the seein : ” —_ ° 
; fined, On the other hand, if the condition is one on the heat 


latter term _ . . . 

Se transfer, then B cannot be found Fortunately, however, this 
n the he pecans i : i ; 
i uncertainty affects neither the shear nor the heat transfer \s 





was pointed out before, no condition far from the plate is of use 


because the effect of this term decays exponentially in the usual 





j rat 
over a flat for large distances In this connection it might be ob 
ilmost 5 i that, for an insulated plate and a constant Pr of unity, 
nding ares } Fin) in Eq. (3) becomes simply 

istant spe- F(n) = Byotton/(pottor 


M2L75 @=75° 
h=37 


~~ 






arbitrar . ° : ‘ 
agg To the extent that this term is nonzero, the temperature distri 

re proniles eh P > . Bae : i 

liffers from the familiar exact integral of the energy equa- 




















bution ¢ 
' tion—namely, 
iijia =i , 1)/2| M2 (1 — u?/u,? re) 
| nt re) 2 4 6 8 10 
2 Phis ambiquity is not related to the slip phenomenon but really x/L 
{ ecurs in the continuum hypersonic analyses (reference 3, for . , ; ; 
{ , Ithough tl uthor eae Gti ak ten a Fic. 1. Comparison between pressure distributions obtained 
mpl Ithoug! 1 g Oo appare au oO otice . “ oe : 
a, ee 2 a Ji saps, 5 isi from characteristics, experiment, and approximate method 
Stewartson? found a similar effect arising in terms of the next order 
' higher than considered in the present study He examined only 
the momentum equation after reducing it to the incompressible ( = ( "os = ( l 
case by the familiar Howarth transformation However, it is seen : ; ; 
j : das ms where ¢ is the two-dimensional pressure coefficient and ( 
here that compressibility effects introduced by the energy equa- 
tion cause the trouble to arise even sooner. It appears that some the conical pressure coefficient for the cowl le In order t 
po] X ; tream condition is required in order to choose the appropriat« satisfy Eqs. (1) the pressure distribution on the cowl surface 
value for B if a heat-transfer boundary condition is used. In any is assumed to be of the form 
} case, there is no effect on the shear and heat transfer although, C,=x2 A Be~# 9 
—_ ‘ for example, the recovery factor cannot be found ’ : ; : 
a Sa ‘ ne where a is a free constant and is the cowl radius at the lip 
The slip effect (which vanishes if a; and c; are set equal to a 5 : 
: ' . Phen, from Eqs. (1 
ro) does not influence either the shear or the heat transfer t 
this order However this conclusion holds only for the flat- Cp, = C a ( cu 3 
} plate case The shear is increased over its continuum valu ere ; : : 
ss 3 : : , : Setting x’ = x/L, A = ( C,.,andrAX = L me obtains 
his effect has nothing to do with slip but is caused by the ex- 
} ternal pressure gradient induced by the boundary layer The Cc, =C A A i 
heat transfer is unchanged from the continuum value. The heat = : _ 
, . . , The constant a may be determined by differentiating Eq. (4 
hickness, { flow due to the sliding friction exactly compensates for the re- ‘ : F 
bata § . : : 2 ; with respect to x’, resulting in 
\/ 72 duced conduction to the wall 
t % ' d C,,/dx = av de 5 
itions i 
; NCES : : : ; : _ - , 
te } REFERENC! If this derivative is evaluated at x’ = 0, the following result ts 
ns have 4} laslen, S. H., Second Approxim ati n to the Laminar ( pressible Bound ee 
H Slip Flow, NACA TN 2818, November, 1952 
it trans- * Martino, R. L., Heat Transfer in Slip Flow, UTIA Report No. 35, Uni a = |d¢ dx Q AA 6 
d Re- versity of Toronto, October, 1955 
it eff I 10, ¥. H., Viscous Flow Along a Flat Plate Moving at High Supersoni 
Spe Journal of the Aeronautical Sciences, Vol. 23 2, pp. 125-13¢ 
lowed | shia f the Aeronautical Science Vol 3 PI 1 14 r = . : = ‘ 
ebruary Yort 
om? i -I19° 
ed limit tewartson. K.. On A mptotic Expansi in the Theory of Bounda M=200 8=10 
} I |. Math. & Phys., Vol. 36, Ne. 3, pp. 173-191, October, 1957 12 T T T T T T T 
t « Dd 
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An Approximate Method for Determining the ff 
0 (4 Wave Drag of Axisymmetric Conical Cowls 06 f | | | 
| I». — METHOD OF 

pe he-issees ee re 04) T +— CHARACTERISTICS 

















} 
. Grumman Aircraft Engineering Corporation / 
02 | // . APPROXIMATE ___ 
METHOD 
r IS ASS h he strear . } _ ic 0 
I IS ASSUMED that the stream tube captured by the cowl 1 O | 2 3 4 5 6 7 a 
re ' vlindrical and that the flow is uniform and steady. There- = 
fore, at the cowl lip, the flow is taken to be two dimensional A=L/r, 
: conical flow will be reached asymptotically as the cowl Fic. 2. Comparison of wave drags obtained by characteristics 





length increases indefinitely Thus and approximate method 
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In order to determine this constant, it is only necessary to de 


termine the derivative |dC,/dx'|,”=9 since \ depends on the 
geometry and A may be found readily This derivative can be 
evaluated by using a simple three-point method-of-characterist ies 
solution The value of the derivative will be a function of the 
mesh size of the characteristic net, but it has been found that it 


does not vary appreciably if small mesh sizes are used 

rherefore, it is seen that the determination of the pressure dis 
tribution on the cowl has been reduced from a lengthy character 
istics solution to one involving only three points. Furthermore, 
the constant a@ may be initially computed for several families of 
cowls and presented in chart form, thereby allowing specific con 
figurations to be calculated rapidly and easily 


sv integrating the known pressure distribution over the cowl 


surface, the wave-drag coefficient may be obtained In evlin 
drical coordinates, the equations of the cowl are 

= T 6x, d = 6d. 7 
Here, 6 tan 6, where @ is the cowl angle For a body of revolu 


tion, the drag coefficient may be written 


. 
Co= (2 ira s 
‘ 
based on maximum frontal area = 7 Phus, substituting 
Eqs. (4) and (7) into Eq. (8 
e/ 
( = (2 c. A éx! 6d. y 
ae () 
Integrating and employing the notation \ = L 
Cn 26/(1 + OA i AC + (1/2)6d2¢ 
A/a) {1 ( ‘) + (46/a?) [e ‘(1 + ad 1} 10 
Phis equation is valid for all values of AX between 0 and , and 


} 


it can be shown to be valid at the two limits Some values ol 
Cy have been plotted in reference (1) for various values of 4, 1/, 
ind 

Phe present method ha been compares to the complete 
method-of-characteristics solution for 1/ 2.00, 0 1O°, and 
\ |, 2, and 5 \lso, the method has been compared with « 
perimental results obtained in) Convair Daimgertiel Dexa 
wind tunne In both cases the greemMent Is gor lig l 
ompares the pressure distributions obt ed using the approxi 
mate method with those from the method of characteristics and 
the experimental values Fig. (2) compares the wave-drag vari 
ition obtained using the present method with that obtained 
from the method of characteristics 
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The Determination of the Shear Center for a 
Special Solid Symmetrical Airfoil 


Alphonse J. Sistino 

Analytical Engineer, Fairchild Engine Division, Fairchild Engine and 
Aircraft Corporation, Deer Park, N.Y 

March 12, 1958 


SYMBOLS 


2 Cartesian coordinates as shown in I l 
FP a force in the xy plane and parallel to the x axi 
shear stress on the xy plane with the force in the x direction 


shear stress on the plane with the force in the y direction 


¢ St. Venant stress function 


T 
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irea Moment ¢ tia abe the ix 
P on's ratic 
| ighses stuby of the deflections of lid airfoil cant 
such as a compressor or turbine blade) the lo f 
shear center plays an important role In gen 
centers for svmmetrical solid airfoils are not readily 
to the difficulty of finding exact analvtical uti 
see below In this paper a special symmetric 
has been chosen such that Eq. (5) will have 
solution his result may serve to check approxi 
for obtaining shear centers of airfoils In Fig t 
foil and its shear center is plotted . 
The shear center is defined by the following equati 
7 ly j> 
d ] / 
Where i JS T T dx 
For a concentrated loaded cantilever beam, t 
equation leads to the following sh stress distributi 
T Oy VO r 21 | 
T ley OV } 
where Isa function su t 
F »] t} 
t the boundar ftheer SEC 
By « dering the compatibility equation, t 
at LIS cer eC 
d*¢ 0. 2 O I P / 
Now if an airfoil section sen whose bou 
given by 
then, bv Eq. (4 
where #1 rbitrary constant to be determine ) 
By substitution of Eq. (6) int 1. (9), Eq. (5) b 
= 
D746 / 0. a7 0 
) 
Now since ¢ must vani t 1 
iq ‘ { 4 Ik 1 [1 l1 1e¢ — 
, i 
where / 4 , & 
Ii stituting t imed ¢ function into Eq. (7 
stants / ind re s ljustes to satisfy Eq. (7 For 
0.3 the f values satisfy Eq. (7 [~ 
/l = | | 
i 
0.0045 (P/J] 
Having now completely defined nd ¢ (x4 T T | 
= 
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Eqs. (2) and (3), and the results are substi TABLE 2 
; wp ] 4) m1 ‘ fio lh 
For h = 0.202, d = 0.061 In Fig. 1 the Incomplete Beta Function Ratios I, (a, 8 
ion is plotted and its shear center noted - > — — ‘ 
REFERENCI 
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Fundamental Solutions for Heat Transfer 
| From Nonisothermal Flat Plates These fundamental solutions involve thi p and ramp di 
tributions defined by Table 1 Phese tutions are also useful 
d and W Reynolds s building blocks either yr ting the heat flux Bt 
Nat 3/ Researct Ti f Inada ttawa, and Stantord it where the wall temperature 7 F. is prescribed or in pre 
/ ; Respective dicting the wall temperature when the heat flux is given. hh 
1958 addition, certain integrals which appear repeatedly in_ the 
problems are tabulated 
ae EARS ther heen Loci, terest in the For problems of specified wall temperature, the heat-transtfet 
. | a ‘ ' temp ture distributi ons rate may be written as 
| the ¢ \ t I f p 
{ fl te { ) ib] . : . 
while for cases of specified heat flux the wall temperaturs 
1 t Vv te rature " 
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Fic. | A wall temperature distribution 


TABLE 3 


Laminar, Turbulent, 


Characteristic SUPr*/sRe, StP/-4Re,°-? 
Constant temperature 0.332 0.0296 
Ramp temperature from leading edge 0.535 0.0336 
Constant heat flux 0.455 0.0309 


Ramp heat flux from eading edge 0.598 0.0342 


for turbulent flow 


» 


Ste = 0.0296 Pr~°-4Re “2 


Here Pr denotes the Prandtl Number, Pr = yuc,/k, and Re 
the local Reynolds Number based on distance from the leading 
edge, Re, = pl'x/p. Fr(x) and F,(x) are the wall-temperature 
and heat-flux functions which describe the nonisothermal in- 
fluence, values of which are given by Table 1 

The solutions in Table 1 are written in terms of ratios of in- 


complete beta function defined by 
Ia, B) = Bia, B B,(a, B 


where the incomplete beta function is defined as 
. 


1 — 
Bia, B) = | o* ‘(1 — ao)? ‘da 
J0 


The complete beta function is related to the gamma function by 
Bila, 8) = I(a)I'(8)/T(a + B) = B,(8, a 

The arguments of the ‘nxcomplete beta function may be shifted 

using the relation 


Bia, B) = Bia, B B,—2(8, a 


The pertinent functions have been determined numerically and 
are presented in Table 2 
For the simplest cases the results in Table 3 may be obtained 
Any problem where the wall temperature or heat flux is 
specified may be handled simply by replacing the given dis- 
tribution by a succession of steps or ramps, and superimposing 
the effects of each This is permissible because the energy 


equation is /inear.’ For example, consider Fig 1 For laminar 


flow the superposition results in 


a"(x) = 0.332 pUc,Pr Re b, <x £4 @ 
q”( *) = 0.332 plc,Pr Re bh + 1.612 mxI (2 3, $/3) > 


Specified heat-flux problems may be treated in a_ similar 


manner 
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Reduction of Ikenberry-Truesdell Equations 


to Burnett Equations for Slip Flow? 


Hsun-Tiao Yang* 
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4 pen PURPOSE of the present note is to establish the « 
equivalence of the second iterates of the Ikenberry-Trues 


Equations! with the Burnett Equations? * for slip floy I 
second-order expression for the stress tensor given by | q. (2.18 
of reference 1 is 
Pj v= — QuE;, T 2u/p we T 2u?/p EE 
$u?/f EA ;u — $u?/3p) Bik 
M/ pO) [3(ue — (p60 
where 
Pe Pij — péb;; = stress tensor 
v coefficient of dvnamic viscosity 
p hydrostatic pressure 
p density 
7] temperature 9 
u gas velocity vector 
E ul Ou}. / OX 
E 1/2 Bij TU = 3) ul 6 OHu,;/O 
D Dt O/ ot + ul 


and a parenthesis around s subscripts indicates a sum over 
s! permutations of the indices, divided by s!. Thus, Eq 


be rewritten in more familiar notation as 


! ) { 
P.\) = —9 Ou = wl Ol 
po a Se 2 M t 
Ox P Dt Ox 
: Bw” Ou, OU : [Ti (~ Ou Ou; Ol ) | 
2 + 2 | : 
p Ox ON, Dp On Ox On On 
tu? Ol; OUK - nh OO O8 
5; +3 yu 
Sp Ox; OA p@ Ox On 
\ccording to the kinetic theory of gases, viscosity depends 
ontemperature. Noting also that 
"id = —2u(Ou;/Ox 
is the first-order Navier-Stokes approximation for the stress 
tensor, one rewrites Eq. (3) as 
pu? D Ou uw? 6 du Dé Ou 
P;9 — Pj) = 2 + 2- 
: , p Dt \ox b » de Dt Ox, 
ms bw? Ou; Olt Pe (~* Ou Ou, Olt ) 
p Ox, On t Ox; OX Ox; OX » 
fu? Ou; OU _ 0° uw? 6 du O68 O86 
0 T » T ks 
op OX; OX p@ Ox;Ox p02 uw d@ OX; O 
Now, all terms in Eq. (5) are of the second order—namely, 
uw’. To be consistent, one introduces the zero-order Euler Equa 
tions to compute (D/Dt) (Ou;/Ox;) and Dé/Dt in the first 
second terms of Eq. (5 From the momentum equation 
Ou/Ot) + uic(Ou;/Ox;,) + (1/p) (Op/dx;) = 0 
we have 
D/Dt\(ou;/Ox;) = — (0/Ox l/p) (Op/Ox;)| — ti 
Ou; /OX;)(OUu;/O ( ° 
S 
From the continuity equation, \ 
+ This work was carried out under Contract AF 18(600)993, sponsored by 
the Office of Scientific Research, Air Research and Development Command h 


and under the supervision of Professor J. M. Burgers 
* Research Associate 
The author wishes to thank Dr. F. S. Sherman for indicating 


minor ( 


errors in the computation 
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Op/Ot) + (0/Ox pu = () ind the equation of state, with R denoting the gas constant 
uati 
1 ons P : t pké, we have 
e energ quation for a monatomic gas, 
i De, Dt — (2/3)0(Ou,/ Ox 7 


Ol + Up(OP/ OX + 2/3) P(OUK/ OX = ({) 
The fourth and fifth terms in Eq. (5) may be combined as 














2(p*/ ft Ou; / OX Ou; / OX t Ol, / OX Ol, / OX _ fu? Sf Ol; / OX OM, OX) )¢ 
$(p*/ ft Ou ¢/ OX Ou; / OX t tu*/OF Ole / OX Ou; / OxX))6 
le exact g 
’-Truesdi 4 (yr? /] Ou,/ OX Ou; / OX 
How l u?/p) |8(Ou,/ Ox Ou, / O2 + (8/3) (Ou; /OX Ou ;/ OX);) — 4(0U;/ OX Ol; OX 
Eq 2.18 
ince O ) () 
Substituting Eqs. (6 8) into Eq. (5) and rearranging terms, we finally have 
y } (: 6 =) w- Ol, OU , we 0 (’ Op ) Ou, Ol : Ou; Ol 
; 3 \2 uw dé] p Ox, OX f Ox; \p Ox Ox; OX Ox, Ox; _] 
_w- (O8 u- OP O8 6 du pw- O68 O86 um” Ou; Oi 
~ () - 3 ~~ & a 
p@ Ox;,OX, ppé Ox Ox u a6 pé- Ox, On f On Ox 
which is exactly the Burnett expression for the stress tensor for Maxwellian molecules given in references 2 and 3 
The second-order expression for the heat flux, /,;, given by Eq. (11.1) of reference 1 is 
2 apd; + (Bap/2p) (ué., + (ap?/2p) (50iE + 36.;u;.; + 86 ,;E + Biu/p) (uE -— 
s(u2/pp) E;;p + L5y/2p8) (uA8 + 12u ip p?/p)E VE + (Spu/it u?/p)E VE + (Op? /ppiEi; (uk 10) 
dA 
ere 
on 15/4)R the slip flow regime than the Burnett equations However, to 
. over t 2/9) Inf ipply the Ikenberry-Truesdell equations to problems in rarefied 
—p - VP 1 } } 
. ()n gas dynamics, the appropriate boundary conditions must first be 
ting that developed 
ht — ap — (15/4 Ru O68 / Ox 1] 
; ; REFERENCES 
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rewrite Eq. (10) as, keeping only yp? terms, Wang Chang, C. S., and Uhlenbeck, G. E., O P} 
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VxB= R,,j 
Phe dimensionless numbers k, &,, and re 


‘a 2) y 
R p P, IX M,o Lo, 5 0 





Here w, is the magnetic permeability, and & is the det 








iecross the shock, which is the only parameter in the ¢ 

field R is called the magnetic Revnolds Number, i: 
with the usual viscous Reynolds Number It is diset 
example, in references 2 ind 5, where it is shown to be ! 
ir less for hypersonic flight conditions Sis a frequently 


is Vet nameless quantity representing the ratio of the mag 


body force to the dynamic pressure.* Unlike &,,, it cont 
field strength and so its magnitude depends on the applic 
For the practically important case of small R,,, Eq. (3 


written V x B = 0), indicating that the currents flowing i 


do not affect the magnetic field, and thus the applie« g 
field is undistorted by the flow Therefore, the B field 
considered given, greatly simplifving the problem Furt 


considerations of axial symmetry combined with Eqs. (2 


3a) show that there is no E field provided the body is 


or an equipotential surface hese two results have 





completely uncoupling Eqs. (1) from Eqs. (2) and (3 i 








, df due toa m field is investig for larg g vhere B, is the app \ t gi f s 
netic Revnolds Number &,, and is found to lead t I sof point consists of o1 me comp 
heat-transfer rate up to 30 per cent However, the practical satisfy T-B 0 it must have the form B ( \ 
interesting case for hypersonic flight is the opposite one of sma constait 
R,,, which is mentioned but not considered in refer r » In Phe pressure 1s eliminated bv taking the cur f Eq 5 
reference 4 a similar problem ts investigated for incompressible stream function identical to thill Is Litroduces 
flow, but oaly shear stress is found. Unfortunately, t results Eq. (ob When B, and thi ( functi ire pu 
of reference 4 are clouded b use tf iuthors compar ist curl of Eq. (Sa) the following differential equat sult 
with a magnetic field to the case without a field but wi the same 

} 2 sin 8 cos OW (d o y 

H d essure dist (tor Phis comparison is not 1 isti 
since for given flight conditions, the field deinit ( ges the oak si ils v 
pressure distribution, as will be shown below Phe primes denote differentiation with respect to 

The purpose of this note is to present the effect of a magueti In the nonmagnetic case S 0, 6) be es 
field on the inviseid flow at the staguation point of a blunt axi linear equation for W for all 6 and can be solved expli 
symmetric body in the practical case of small magnetic Reynolds In the present case it becomes an ordinary nonlinear equat 
Number In view of the conclusions of reference 2, that the main the stagnation point, 9 () 


effect of the field is to change the inviscid rather than the viscous 


flow, this will make possible a prediction of the change in heat 





transfer and pressure drag caused by the field The method used 


v Lighthill,® where a 


to analyze the flow is that introduced | 
spherical shock shape and incompressible flow are assumed 0.94 
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The blunt body is shown in Fig. 1 in a stream of velocity / 
and density p Near the stagnation point it is taken as a sphere 
of radius ”,*, producing a spherical shock of radius 7,*, behind 0.7 + 
Which there is an incompressible flow of density p and constant 


conductivity o. The applied magnetic field B,* is taken to be in 06 





the radial direction only and spherical coordinates ¢* and @ are 


used 

The governing magnetohyvdrodyvnamic equations relate the 05" 
magnetic field B*, velocity g*, pressure p*, electric-current density 
J*, electric field E*, and vorticity Q* = V x q* They may be Ss = = s | 4 
found, for example, in references 2 or 4. If the equations are 


made dimensionless, using a length Z, velocity l'.., density p 





field B,*, and conductivity o, we find the following steady-state 
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equations for the dimensionless quantities (without asterisks 
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» rs: _ 2 > 
Vikp + q?/2 qx2 = SkjxB la S_p=o(B,)" /RQUS 
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g Lig i 
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I ) Phe re 1 A 
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siles rhe pressure increase is somewhat larger, being pr entrifugal fore ( 
il to (1 S,/4 Phis factor is also shown in Fig. 2. A wv the r pressure gr If t 
stunate of the pressure-drag change, obtained by integrat ced t ger 1 fugal f ru 
Eq. (8) for a sphere up to the angle where vanishes, shows tic re es \ \ 
( ig increases by a factor |(1 Sp/4) 1 \ For f d tangent \ ( f f 
le case quoted above, where & = 1/15, this factor is 1\ 
igain the pressure-drag increase is quite small, thoug ‘ ages thon a A i 
( ulation is admittedly very crude since the present theory 
t extend far from the stagnation point Allin all, it would 


pear that the sort of 


geometry considered here is suitable 
ining significant effects near the stagnation points of 
Phe situati Is, Of course, more favorable for larger 
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the angular momentum of a fluid particle is preserved so that 
roy = 7, Uy, and the centrifugal force is now p;[(71? u,2)/r23 If 
this force is larger than the prevailing pressure gradient 
p2(o2/re), or if Py 42 uy;2 > ps» U2? ro2, the fluid element will con- 
tinue its motion outward and the motion is unstable. On the 
other hand, if p; 71? 1? p2 U2? ro*, the motion is dynamically 
stable 

In von Karman’s original analysis the density of the fluid is 
uniform, so that pi = p2 = p:, and the famous Rayleigh criterion 


is obtained—i.e., 


dk?/dr < 0 unstable 
> () stable 


where k? = (ru)? is the square of the angular momentum. In 
general, the new density of the fluid particle p; is almost cer- 
tainly different from p;, but the strongest destabilizing effect 
Then the 


is obtained if we take p: = pi, because p, p» here 


new criterion reads 
d/dr(p u2r?) <0 unstable 
> 0 stable 


or 


unstable! 


p ur242[(O0u/Or) + (u/r)| +(u/p)(Op/dr)} 
stable | 


In the case of the laminar boundary layer over a highly cooled, 
blunt, convex body, the radial coordinate is r = R + y where 
y is the distance normal to the surface and 6/R < 1. (See 
Fig. 1.) Now the density in the boundary layer is given ap- 
proximately by p./p = h/h,., where h is the enthalpy, or 


(u/p)(Op/Or) = —(u/h)(Oh/odr 


But when p./pw < 1, a good approximation to the enthalpy 
distribution is given by the simple relation h >= h, + (hk —- 
h,)(u/u,.) where the subscript w refers to surface conditions, and 
the subscript e refers to conditions at the outer edge of the 


boundary layer. Therefore, 


(u/h)(Oh/Or) = (u/h\(Oh/oy) = 
\(te — Iw)(u/ue)/[Aw + (he — hw)(u/u)] { (Ou/dy) 


and the quantity in the braces in Eq. (1) becomes 


ou J, (he — hwu/ue) 2u 


oy t- Iw + (he — hw )(u/ue)§ + y 
Now the quantity 
(he — hw )(t/tte)/ [hw + (he — hw) t/t 


is never larger than unity, so long as h, < h,, and the quantity 
2u/(R + y) is always positive (but negligibly small Also 
ou/dy >0. Weconclude therefore that the criterion for dynamic 
stability given by Eq. (1) is always satisfied in this case 

It must be emphasized that this analysis is not meant to 
apply in the immediate vicinity of the forward stagnation point 
where the streamline curvature of the inviscid flow outside the 
boundary layer is concave, and the dynamic instability may be 


significant 
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t This criterion is identical with the result obtained by Synge? for liquids 
by applying the Sturm-Liouville oscillation theorem to the differential 
equation for small disturbances in a laminar motion In his case dp/ dt 0, 
and 7 is in fact equal to p An analysis for a gas with arbitrary Mach 
Number is now being made but is not expected to change these results very 
much. This result can also be shown to be identical with Eq. (35) of the 


recent paper by M. Lessen.3 
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An Evaluation of the Criticisms by Topping of 
Marin’s Strength Hypothesis 


John J. Coleman 

Mathematical Analysis and Operations Research Group 
Electro Dynamics Staff, Transport Division, Boeing Airplane 
Company, Renton, Wash. 

March 7, 1958 


gap iposeeoi MARIN STATED that his proposed theory of failure 
of materials under combined stress was a generalization of 
the octahedral shear stress theory” 


° , ) , , 
Gi oe” SE Ha? oe — Si)” (1) 


where 01, 92, 03 are principal stresses and 7» is the octahedral 
shear stress 
Marin’s relation states 


(o, — a)? + (o2 — 6)? + (¢ 


(o2 — b)(a3 — ¢ 


where og,” is a constant of the material, and a, b, c, and gq are also 
constants 

In his note,! Topping tried to show, by demonstrating its 
shortcomings, that Marin’s proposal is untenable for various 
values of g. Unfortunately, all present theories of failure for 
combined stress have shortcomings of one kind or another. 

Topping’s insistence that Marin’s Eq. (2) fit all the restrictions 
imposed from known behavior of materials under all conditions 
of failure ‘seems overcritical. However, the imposed restrictions 
mentioned by Topping do restrict the value of g to less than —1, 
In general, Marin’s formula then satisfies the conditions for 
unequal tensile and compressive strengths, anisotropy, and the 
features that the limiting failure surfaces must be open at the 
compression end and closed at the tension end of the hydrostatic 
stress axis. The last set of features also include the fact that 
failure is affected by hydrostatic pressure 

Strictly speaking, if all points on the yield surface represented 
by Marin’s Equation for gq < —1 represent the condition of 
vielding, it appears to contradict the experimental facts because 
the state of hydrostatic tension of three equal tensile stresses, 
under which a material cannot yield, would correspond to such 
a point. However, a second surface, distinct from the yielding 
surface, must exist—-namely, the surface of rupture. One can 
assume that the yield surface encloses the rupture surface 
beyond a common intersection curve in the positive quadrant. 
A geometrical situation such as this would account for rupture 
occurring before vielding in the case of triaxial tension. Conse- 
quently, the hyperboloid of one sheet, defined by g < —1, may 
be applicable after all. 

In conclusion, one must say that Topping performed a service 
for Professor Marin in determining the restrictions on g. How- 
ever, Topping’s arguments did not make Eq. (2) untenable for 
the case g < —1. In fact, this case demands the existence of 
two limiting surfaces for triaxial tension. Similar mathematical 
expressions might be used to define both. The seriousness of 
any shortcomings should be judged in the light of the usefulness 
of the proposed theory in providing adequate design criterion. 
Furthermore, the proposed theory should be compared with test 


results before any final judgment can be made 
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